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Three  large  computer  codes  are  included  in  the  Surface  Wave  Analysis 
Package.  The  first,  TELVEL,  obtains  phase  and  group  velocities  and  spectral 
amplitudes  from  a  seismogram  by  using  narrow-band  filtering  and  phase-matched 
filtering.  The  second,  INVERT,  inverts  the  phase  and  group  velocities  for 
earth  structure,  and  then  inverts  for  moment  and  Q  structure,  by  comparing 
synthetic  and  observed  spectral  amplitudes.  The  third,  SYNSRF,  generates 
synthetic  seismograms,  dispersion  curves  and  eigenfunctions.  These  three 
codes  are  used  together  to  obtain  a  Green's  function  for  the  source-station 
travel  path. 

I 

For  each  SRO  station  travel  path  we  obtain  phase  and  group  velocities 
and  an  average  path  structure.  In  addition,  the  narrow-band  filtered 
seismogram  provides  an  indication  of  station  quality  by  revealing  stations 
!  with  severe  multipathing  or  other  problems.  Stations  KONO  and  SHIO  prove  to 
be  excellent  stations  with  very  clean  group  velocity  curves.  Stations  MAJO, 
GRFO,  ANTO,  and  CHTO  show  clear  evidence  of  multipathing,  but  are  still 
useable.  Station  KAAO  shows  severe  multipathing  resulting  in  a  bifurcated 
group  velocity  curve. 

This  report  is  intended  to  provide  a  users  manual  for  the  Surface  Wave 
Analysis  Package.  The  appendices  contain  file  formats,  sign  conventions,  an 
explanation  of  some  of  the  algorithms  used  in  the  programs,  and  definitions  of 
output  quantities,  as  well  as  a  detailed  description  of  program  TELVEL. 
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I.  SURFACE  WAVE  ANALYSIS  PACKAGE 


This  report  describes  the  S-CUBED  Surface  Wave  Analysis 
Package  which  is  designed  to  obtain  path  corrections  from  observed 
seismograms.  There  are  three  main  programs  in  the  package.  The 
first  —  TELVEL,  uses  the  method  of  phase-matched  filters  (Herrin 
and  Goforth,  1977)  to  recover  surface  wave  phase  and  group 
velocities.  The  second  —  INVERT,  inverts  the  phase  and  group 
velocities  to  obtain  the  earth  structure  for  the  path,  and  obtains 
an  approximate  Q  structure  and  moment  by  comparing  synthetic  and 
observed  spectral  amplitudes.  The  third  —  SYNSRF,  uses  a  variety 
of  methods  (Takeuchi  and  Saito  (1972),  Schwab  and  Knopoff  (1970), 
Harkrider  (1964,  1970))  modified  for  high  frequency  stability,  to 
synthesize  surface  waves. 

All  programs  are  interactive  and  self-contained.  Every  effort 
has  been  made  to  make  them  as  easy  to  use  as  possible.  In  this 
report  we  show  tne  uses  of  the  programs  for  making  path  corrections 
by  going  through  an  example  in  detail  using  a  synthetic  seismogram, 
and  by  applying  the  method  to  Shagan  River  —  SRO  station  travel 
patns. 


1.1  WHAT  IS  A  PATH  CORRECTION? 

A  patn  correction  is  a  Green's  function  for  a  given  source 
region  and  source-to-receiver  path.  Once  the  Green's  function  is 
known,  it  may  be  used  to  generate  synthetic  surface  waves  from  an 
arbitrary  source,  as  observed  at  the  receiver  point,  or  the  inverse 
problem  may  be  solved:  an  observed  seismogram  may  be  used  to 
estimate  the  strength  and  type  of  the  source,  usually  in  the  form  of 
a  moment  tensor. 

The  basic  procedure  used  to  make  a  path  correction  was 
outlined  in  an  earlier  report  (Wang,  et  aj^. ,  1981).  Since  then  the 
procedure  has  been  improved  by  allowing  simultaneous  estimation  of  Q 
and  moment,  and  by  making  the  codes  interactive  and  compatible. 


1.2  HOW  TO  MAKE  A  PATH  CORRECTION 

To  make  a  path  correction,  simply  take  a  seismogram,  use 

program  TELVEL  to  obtain  the  Rayleigh  wave  phase  and  group 
velocities,  then  use  program  INVERT  to  find  the  average 

source-to-receiver  shear  velocity  and  Q  structure,  and  finally  use 
program  SYNSRF  to  generate  eigenfunctions  for  this  structure,  and  to 
compute  a  synthetic  seismogram  which  resembles  the  original 

seismogram.  The  scalar  moment  of  the  explosion  is  found  together 

with  the  estimate  of  Q,  and  the  eigenfunctions  can  be  used  to 
perform  a  moment  tensor  inversion  for  any  source  in  the  same  area  as 
the  original  explosion. 

In  practice,  finding  a  path  correction  is  more  than  a  routine 
operation,  and  there  are  a  number  of  ways  to  err  during  the  data 
processing.  The  purpose  of  this  section  is  to  go  through  the 
procedure  in  some  detail  and  to  identify  the  traps  which  are  lurking 
to  confuse  the  user.  Fortunately,  most  of  the  errors  which  can  be 
made  will  show  up  somewhere  during  the  processing.  A  bad  set  of 
phase  and/or  group  velocities  may  not  produce  a  reasonable  earth 
model  when  inverted.  An  inaccurate  inversion  for  shear  velocity  may 
make  it  impossible  to  determine  a  reasonable  Q  structure.  If  the 
seismogram  ootained  at  the  end  of  the  procedure  is  a  good  match  to 
the  original  (except  for  noise  or  multipathing  effects,  etc.),  then 
chances  are  very  good  that  the  path  correction  has  been  accurately 
found. 

The  fii'st  step  in  making  a  path  correction  is  to  collect  all 
of  the  information  pertaining  to  the  seismogram.  The  following 
items  are  required: 

1.  A  seismogram  sampled  at  evenly  spaced  points  in  REALIO 
format  (see  Appendix  1),  preferably  demeaned,  detrended 
and  tapered 

2.  The  number  of  points  in  the  seismogram 

3.  The  sampling  interval  for  the  seismogram 
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4.  The  time  delay  between  the  source  time  and  the  start  of 
the  seismogram 

5.  The  distance  from  the  source  to  the  receiver 

6.  The  instrument  response  expressed  as  a  ratio  of 
polynomials  also  in  REALIO  format 

It  is  very  important  to  have  an  accurate  estimate  of  source  to 
receiver  distance,  and  the  seismogram  start  time.  It  is  well  worth 
double  checking  these  numbers  since  everything  that  follows  will  be 
wrong  if  these  numbers  are  incorrect.  Of  course,  errors  of  a  few 
seconds  or  a  few  kilometers  may  not  result  in  major  errors  in  the 
final  results,  but  it  is  the  abundant  typos  and  minute  errors  that 
can  be  disastrous.  The  best  way  to  be  sure  of  timing  is  to  compare 
several  seismograms  for  the  same  path  to  be  sure  the  times  are  all 
consistent. 

Several  small  utility  codes  are  available  which  allow  the 
conversion  of  any  seismograms  to  REALIO  format  as  well  as 
detrending,  demeaning,  etc.,  and  which  put  instrument  responses  in 
the  proper  format  (see  Appendix  2).  Instrument  responses  must  be 
known  as  polynomials;  no  provision  is  currently  made  for  tabulated 
instrument  responses.  One  additional  quantity  is  required  by  TELVEL 
—  the  initial  phase  of  the  source.  For  an  explosion  (vertical 
component,  displacement,  positive  up,  Rayleigh  wave)  the  initial 
phase  is  -  3tt/4.  In  order  to  illustrate  the  procedure  we  apply  the 
entire  surface  wave  analysis  package  to  a  synthetic  seismogram. 
This  is  a  useful  exercise  to  go  through  before  processing  data  for  a 
new  area,  as  the  results  help  to  identify  certain  problems  in 
advance,  and  may  help  to  identify  the  correct  Zn  branch  in  the  phase 
velocity  analysis.  If  an  approximate  structure  is  available  for  an 
area,  it  should  be  used  to  make  a  similar  test. 

A  synthetic  seismogram  was  constructed  (with  program  SYNSRF) 
using  the  structure  listed  and  shown  in  Figure  1.  The  synthetic  was 
computed  using  an  SRO  long  period  instrument  at  a  distance  of  3000 
kilometers  and  an  explosion  source  with  a  (the  long  period  limit 
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THICKNESS 

(m) 

ALPHA 

(m/sec) 

SETA 

(nj/sec) 

RHO 

(Kg/m^) 

Q 

1100.0 

5000.0 

2700.0 

2100.0 

150.0 

25000.0 

5900.0 

3200.0 

2500.0 

250.0 

19000.0 

6800.0 

3700.0 

2800.0 

400.0 

20000.0 

8100.0 

4500.0 

3300.0 

600.0 

15000.0 

8200.0 

4500.0 

3300.0 

600.0 

20000.0 

8000.0 

4400.0 

3300.0 

100.0 

100000.0 

7800.0 

4200.0 

3200.0 

90.0 

20000.0 

8000.0 

4400.0 

3300.0 

100.0 

80000.0 

3500.0 

4700.0 

3500.0 

200.0 

Figure  1.  East  Kazakh  shear  velocity  structure  used  for  synthetic 
seisinograin. 
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of  the  reduced  velocity  potential)  of  one  cubic  meter.  The 
synthetic  seismogram  and  t»  ue  phase  and  group  velocity  dispersion 
curves  are  shown  in  Figure  1.  Figure  3  shows  a  form  summarizing  tne 
inputs  and  outputs  to  the  programs.  It  is  useful  to  have  a  similar 
form  available  before  starting  the  procedure. 


1.2.1  Recovery  of  Phase  and  Group  Velocities  and  Spectral  Amplitudes 

With  this  information  in  hand,  the  next  step  is  to  run 
TELVEL.  The  first  phase  of  TELVEL  is  to  apply  a  set  of  narrow  band 
filters  to  the  seismogram  to  obtain  an  approximate  set  of  group 
velocities.  The  choice  of  filter  frequencies  and  filter  widths 
(specified  in  terms  of  tne  filter  Q)  is  up  to  the  user.  From 
experience,  a  Q  of  10  provides  a  good  filter  width  for  most  data. 
For  long  period  data,  sucn  as  SRO  data,  20-40  frequencies  from  .01 
to  .1  Hz  usually  produces  a  good  group  velocity  curve  over  the 
reliable  frequency  range  of  data  (see  Figure  4).  Occasionally,  the 
program  will  encounter  difficulties  in  finding  a  good  group  velocity 
curve.  If  it  does,  try  going  back  and  using  more  or  fewer 
frequencies  or  a  lower  or  higher  Q.  Sometimes  the  program  will  find 
a  good  set  of  group  velocities  as  evidenced  by  the  peaks  on  the 
plot,  but  will  not  recognize  it  (no  curve  will  be  drawn  through 
them).  If  so,  they  can  be  entered  by  hand  using  the  edit  command. 
The  narrow  band  filter  is  an  important  step  in  finding  the  group  and 
phase  velocities.  Not  only  does  it  provide  an  initial  group 
velocity  estimate,  but  it  may  show  up  features  such  as  bifurcated 
group  velocity  curves  which  may  render  the  phase  matched  filter 
unstable. 

The  purpose  of  TELVEL  is  to  improve  the  initial  estimate  of 
group  velocity  by  phase  matched  filtering  while  eliminating 
interfering  multipath  arrivals,  higher  modes  and  noise,  and  to  find 
the  phase  velocities.  The  phase  matched  filter  is  found  by 
integrating  the  group  delay  found  initially  by  narrow  band 
filtering.  When  the  filter  is  applied  to  the  original  seismogram, 
the  true  surface  wave  is  compressed  into  a  narrow  time  window.  The 


Figure  3 


Source  Location:  East  Kazakh 

Receiver  Location:  Synthetic  »  3000  km 

Source  Time:  0  sec 

Receiver  Time:  600sec 

Delay  Time:  600  sec 

Distance:  3000  km 

Seismogram  File;  RSYN  *  EKSRO 
Number  of  Points:  901 
aT:  1.0  sec 
Instrument;  SRO-LP 
Initial  Phase  =  -  0.75 

Narrow  Band  Filters  Used 
NF,  FMIN,  FMAX,  Q  =  30,  .01,  .1,  10 

TELVEL  Output  File(s):  TEL  *  EKSRO 

Eigenfunction  File  for  Source  Region  EF  *  EKAZ2 
Instrument  Gain  (conversion  to  meters)  1. 
Estimated  Moment  =  6.4  x  10^^  /Gain  =  6.4  x  10^^ 
Estimated  »  -980  /Gain  =  .980 

INVERT  Output  File(s):  INV  *  EKSRO 

SYNSRF 

Frequencies  Computed:  100,  .002,  .2 
Output  Files:  EF  *  EKSRO 
Estimated  M^:  0.04 

Comnents:  Test  on  Synthetic  Seismogram 

PATH  CORRECTION  INFORMATION 

Form  summarizing  inputs  and  outputs  for  Surface  Wave 
Analysis  Package. 
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Figure  4.  Group  velocity  curve  obtained  by  narrow  band  filtering  the  synthetic  seismogram  in  Figure 


resulting  phase  matched  filter  output  is  plotted  by  TELVEL  to  allow 
windowing  of  the  time  series  (see  Figure  5).  It  is  important  to 
choose  a  time  window  large  enough  to  contain  most  of  the  energy  in 
the  main  arrival,  but  small  enough  to  eliminate  most  of  the  noise 
and  multiple  arrivals  from  the  filter.  The  best  way  to  choose  the 
time  limits  is  to  look  at  the  phase  matched  filter  output  and  then 
back  up  to  apply  the  appropriate  time  window.  For  long  period  data 
a  time  window  of  50-100  seconds  on  each  side  of  the  main  pulse  is 
usually  appropriate. 

After  the  phase  matched  filter  is  applied  to  the  seismogram, 
the  filter  phase  is  unwrapped  and  used  to  estimate  the  phase 
velocities  of  the  Rayleigh  wave.  Since  any  multiple  of  2it  may  be 
added  to  the  phase,  there  is  an  uncertainty  introduced  into  the 
phase  velocity.  If  ^  is  the  phase  of  the  surface  wave  («5  is  always 
<  0)  then  the  phase  velocity  is  given  by 

-  2irfr 

i  +  2irn 

where  r  is  the  distance,  is  the  initial  phase,  f  is  the 
frequency  and  n  is  an  undetermined  integer.  The  choice  of  the 
correct  value  of  n  is  often  difficult.  In  order  to  help  choose  the 
correct  value  of  n,  TELVEL  plots  the  phase  velocity  for  a  user 
specified  value  of  n  together  with  values  of  n  ±  1  (see  Figure  6). 
The  following  facts  help  in  the  choice  of  n. 

1.  The  true  phase  velocity  approaches  a  constant  in  the  low 
frequency  limit.  If  the  phase  velocity  is  well 
determined,  values  of  n  which  are  too  small  cause  the 
phase  velocity  to  decrease  rapidly  at  low  frequencies. 
While  values  of  n  which  are  too  large  cause  the  phase 
velocity  to  increase  rapidly  at  low  frequencies.  In 
practice,  the  correct  phase  velocity  is  usually  the  first 
curve  to  turn  upwards  at  low  frequencies,  although  it  is 
sometimes  the  first  to  turn  downwards,  and  there  is 
sometimes  no  clear  choice. 
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Figure  5.  Phase  matched  filter  output  for  synthetic  seismogram.  Energy  is  compressed  to  within 
±  100  seconds  of  primary  arrival  time  based  on  initial  group  velocity  curve. 


2.  The  true  phase  velocity  is  not  very  different  from  the 
group  velocity  at  low  frequencies,  and  is  always  greater 
than  the  group  velocity  at  the  lowest  accurately 
determined  frequency. 

3.  The  phase  velocity  is  approximately  equal  to  4  km/sec  in 
most  regions  of  the  world  at  frequencies  between  .01  and 
.02  Hz. 

If  the  choice  of  number  of  2it's  to  add  to  the  phase  is  not 
clear,  save  output  files  for  more  than  one  value  and  see  which  one 
leads  to  a  reasonable  inversion  model.  Also,  processing  of  another 
seismogram  for  the  same  path  may  lead  to  a  clearer  choice  of  phase 
velocity. 

An  iterative  procedure  is  used  to  refine  the  group  velocity 
estimate.  Usually  two  or  three  iterations  will  produce  stable  group 
velocity  estimates  which  do  not  change  from  one  iteration  to  the 
next  (see  Figure  7).  It  is  a  good  idea  to  save  the  phase  velocity 
plot,  and  group  velocity  plot  on  the  final  iteration  for  future 
reference.  The  choice  of  2if's  to  be  added  to  the  phase  may  be 
easier  after  the  group  velocity  has  converged. 

The  final  decision  to  be  made  is  the  choice  of  frequencies  to 
output.  These  do  not  need  to  be  the  same  as  the  initial  narrow  band 
filter  frequencies.  The  range  of  frequencies  should  be  large  enough 
to  include  all  data  points  that  are  reliable,  but  small  enough  to 
exclude  unreliable  data  points.  Two  figures  obtained  in  the 
processing  help  to  identify  the  appropriate  range.  Points  on  the 
initial  narrow  band  filtered  group  velocity  plot  with  amplitudes 
less  than  10  to  20  per  cent  (indicated  by  1  and  2  on  the  plot)  of 
the  maximum  amplitude  are  not  likely  to  be  reliable.  On  the  final 
group  velocity  plot,  portions  of  the  curve  at  the  ends,  which  have 
not  stabilized,  are  not  reliable.  Portions  of  the  group  velocity 
curve  which  appear  odd  such  as  those  which  have  a  rapid  decrease  in 
velocity  or  an  unusually  high  velocity  at  low  frequency  or  rapid 
oscillations  at  high  frequency  should  not  be  used.  Periods  longer 
than  the  time  window  used  are  also  unreliable.  It  is  important  to 
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Figure  7.  Final  group  velocities  obtained  for  synthetic  seismogram  by  iterative  application  of 
phase  matched  filter. 


obtain  information  at  low  frequencies  since  these  frequencies 
determine  the  deep  structure  of  the  earth.  For  most  SRO  stations  at 
distances  of  2000  through  5000  kilometers,  phase  and  group 
velocities  may  be  reliably  obtained  from  about  .015  to  .08  Hz.  A 
convenient  frequency  spacing  is  .005  Hz.  Phase  and  group  velocities 
can  almost  always  be  obtained  down  to  a  minimum  frequency  of  .02  Hz. 

The  output  file  from  TELVEL  contains  the  estimated  phase  and 
group  velocities  and  the  instrument  corrected  spectral  amplitudes  of 
the  seismogram.  The  corrected  amplitudes  are  an  improvement  over 
tne  Fourier  amplitudes,  since  they  are  obtained  from  the  compressed 
seismogram  after  the  removal  of  noise  and  interfering  phases. 

The  TELVEL  output  file  also  contains  a  space  for  the  standard 
deviation  of  the  phase  and  group  velocity.  TELVEL  does  not  estimate 
errors  in  the  velocities,  so  it  simply  outputs  default  values  of 
zero.  If  more  than  one  seismogram  for  the  same  path  are  processed, 
they  can  be  averaged  to  get  an  improved  estimate  of  phase  and  group 
velocities  and  to  estimate  uncertainties  in  the  data.  A  small 
utility  program  (AVEDAT)  is  available  which  constructs  a  combined 
file  in  the  format  of  a  TELVEL  file. 

1.2.2  Inversion  for  Earth  Structure 

The  TELVEL  output  file  is  in  the  proper  format  to  be  used  by 
the  second  major  code  —  INVERT,  wnich  inverts  the  phase  and  group 
velocities  for  the  average  shear  velocity  structure  along  the  path, 
and  which  uses  the  amplitudes  to  estimate  the  Q  structure  of  the 
path. 

Program  INVERT  has  been  designed  to  require  a  minimum  of 
operator  input.  A  starting  model  is  automatically  generated  from 
the  input  data.  The  final  model  is  independent  of  the  starting 
model,  but  it  is  important  to  choose  a  suitable  range  of  depths 
appropriate  for  the  resolution  capability  of  the  data. 

Only  the  TELVEL  output  file  is  required  for  shear  velocity 
inversion.  One  other  file  is  needed  for  Q  inversion.  Assuming  that 
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the  structure  of  the  source  region  is  approximately  known 
independently  from  previous  investigations,  an  eigenfunction  file 
(generated  for  the  source  region  structure  by  program  SYNSRF)  will 
be  input  at  this  stage. 

INVERT  is  divided  into  three  main  sections  —  model 
generation,  inversion,  and  analysis.  Most  of  the  time  the  only 
command  necessary  in  the  model  generation  section  is  RUN  which 
causes  model  generation  to  performed.  In  the  model  generation 
section,  the  user  can  set  a  value  of  Poisson's  ratio  and  the 
coefficients  of  Birch's  law  if  desired.  The  inversion  code  only 
inverts  for  shear  velocity  while  compressional  velocity  and  density 
are  constrained  by  these  constants.  The  default  values  are 
appropriate  for  most  of  the  earth  and  rarely  need  to  be  changed. 
Two  features  which  are  sometimes  needed  are  the  ability  to  set  a 
discontinuity  at  a  particular  depth,  and  the  ability  to  set  the 
number  of  layers  used  in  the  inversion.  The  default  number  of 
layers  is  20.  A  smaller  number  of  layers  will  result  in  a  faster, 
but  less  accurate  inversion. 

In  many  cases  a  discontinuity  may  be  required  to  match  the 
data,  but  its  depth  is  not  usually  known  in  advance.  The  best  plan 
is  to  run  the  inversion  section  and  insert  a  discontinuity  where  the 
resulting  model  shows  a  strong  velocity  gradient. 

The  first  step  in  the  inversion  is  to  constrain  the  smoothness 
of  the  model  defined  by  the  number  of  degrees  of  freedom  (command 
OF)  allowed.  A  small  value  of  OF  produces  a  smooth  model  while  a 
higner  value  produces  a  better  data  fit.  A  good  value  to  start  with 
is  a  OF  of  4. 

The  command  INVERT  causes  the  first  iteration  of  the  inversion 
to  be  performed.  This  command  causes  partial  derivatives  3c/3B  and 
3U/3  9  to  be  calculated  from  the  initial  (or  current)  model.  It 
causes  a  set  of  matrices  to  be  assembled  from  the  partial 
derivatives,  and  performs  a  full  singular  value  decomposition  of  the 
matrices.  Finally  it  computes  the  exact  nonlinear  phase  and  group 
velocities  for  the  estimated  model  with  the  OF  specified  initially. 
All  of  this  requires  quite  a  lot  of  computation  and  will  take  some 
computer  time. 
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After  the  inversion  has  been  performed,  the  command  PLOT  will 
plot  the  shear  velocity  structure  while  the  command  PLOT  DATA  will 
plot  the  linearized  data  fit. 

The  command  OF  has  a  new  meaning  after  an  inversion  has  been 
performed.  A  major  advantage  of  the  matrix  decomposition  is  that 
new  models  for  different  values  of  OF  may  be  generated  without 
recomputing  the  matrices.  Thus  if  the  command  OF  4  was  given  before 
inversion,  the  command  OF  5  after  inversion  instantly  produces  the 
(linearized)  model  for  a  OF  of  5.  The  commands  PLOT  and  PLOT  DATA 
again  produce  the  model  and  data  fit  for  the  new  OF  (see  Figure  8). 
This  makes  it  very  easy  to  pick  the  best  value  of  OF.  A  second 
iteration  can  then  be  performed  using  this  model  as  the  new  starting 
model.  Usually  three  iterations  with  the  same  OF  are  sufficient  for 
convergence. 

An  obvious  question  is  how  to  choose  the  best  model.  The 
answer  is  to  use  a  OF  high  enough  to  get  a  good  data  fit,  but  small 
enough  to  prevent  oscillations  of  the  shear  velocity  model,  and  to 
yield  a  smooth  fit  to  the  data.  Usually  a  OF  of  about  5  gives  a 
good  model.  With  good  quality  data  such  as  the  SRO  data  a  OF  of  6 
may  be  used,  while  less  accurate  data  may  only  allow  a  OF  of  4  (OF 
need  not  be  an  integer). 

One  otner  choice  must  often  be  made  during  inversion.  It  may 
be  necessary  to  include  a  discontinuity  in  the  model.  The  most 
conmon  place  a  discontinuity  is  needed  is  at  the  crust-mantle 
boundary.  The  proper  location  for  the  discontinuity  may  be 
determined  by  the  midpoint  of  a  sharp  gradient  which  appears  in  the 
shear  velocity  model.  Within  the  inversion  section  a  discontinuity 
may  be  added  by  specifying  the  number  of  the  layer  above  the  desired 
discontinuity,  or  it  may  be  specified  by  depth  by  returning  to  the 
model  generation  section  (this  is  unnecessary  unless  the  user 
requests  a  specific  depth  which  does  not  coincide  with  a  layer 
boundary).  If  higher  frequency  data  are  available,  as  it  may  be  for 
WWSSN  instruments  at  distances  less  than  1000  kilometers,  it  may  be 
necessary  to  include  a  shallow  discontinuity  in  the  upper  few 
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Figure  8.  Shear  velocity  structure  obtained  by  inverting  dispersion 
curves  produced  by  TELVEL  from  the  synthetic  seismogram, 
together  with  linearized  data  fit.  The  strong  gradient  from 
20  to  60  km  indicates  a  possible  discontinuity.  This  model 
is  shown  for  the  second  iteration  using  a  OF  of  5.0. 
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l<i lometers.  It  may  be  difficult  to  define  the  proper  depth,  since 
the  pnase  and  group  velocities  may  lose  accuracy  at  the  high 
frequencies  neeoed  to  resolve  shallow  depths.  When  the  presence  of 
tnis  layer  is  indicated,  it  is  often  essential.  The  NTS-Tucson 
seismograms,  for  example,  absolutely  cannot  be  matched  without  a 
Shallow  low  velocity,  low  Q  layer  in  the  upper  two  kilometers 
(3ache,  T.  C.,  W.  L.  Rodi,  and  0.  G.  Harkrider  (1978)). 

After  the  last  iteration  has  been  performed,  give  the  command 
END  to  finalize  the  inversion  and  go  to  the  analysis  section.  Here 
you  can  make  final  plots  of  the  model,  and  the  true  nonlinear  data 
fit  (Figure  9).  You  can  also  print  the  data  fit,  variances  and 
spreads  for  the  final  model. 

One  more  calculation  is  required  to  obtain  the  complete  earth 
model  —  the  Q  structure  needs  to  be  found.  The  program  does  this 
by  computing  a  synthetic  seismogram  (actually  an  unattenuated 
synthetic  spectrum)  at  the  observed  frequencies  and  taking  the  ratio 
with  the  observed  spectrum.  The  logarithm  of  the  ratio  is  equal  to 
a  constant  (the  moment)  plus  a  frequency  dependent  attenuation 
coefficient.  In  effect  the  level  of  the  spectral  ratio  gives  the 
moment,  while  the  shape  of  the  spectral  ratio  as  a  function  of 
frequency  gives  the  attenuation  coefficients  which  can  be  inverted 
for  the  Q  structure.  As  mentioned  before,  if  an  eigenfunction  file 
is  available  for  the  source  region,  the  program  will  generate  the 
synthetic  spectrum  using  these  eigenfunctions  for  the  source  region 
and  using  the  inverted  structure  for  the  path.  The  program  uses  the 
spectral  ratio  and  inverts  for  log  moment  and  s/Q  simultaneously. 

The  spectral  amplitudes  are  never  known  as  accurately  as  the 
phase  and  group  velocities,  so  a  lower  DF  must  be  used  for  Q 
inversion.  The  minimum  possible  DF  is  2.  A  DF  of  2  will  produce  Q 
proportional  to  3  in  all  layers.  An  important  function  of  the  Q 
inversion  is  in  fact  to  smooth  the  spectral  ratios  and  the 
attenuation  coefficients.  It  is  usually  necessary  to  use  a  OF 
between  2  and  3.  Higher  DF's  result  in  unrealistic  Q  structures. 
Discontinuities  are  not  allowed  for  Q  inversion,  and  are  not  used 
even  if  they  were  used  in  the  velocity  inversion. 
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Figure  9 


Final  model  obtained  using  a  OF  of  6.0  and  a  single  discontinuity 
(compare  with  Figure  1)  together  with  data  fit.  The  intermediate 
layer  from  30  to  50  km  has  been  split  by  the  single  discontinuity. 
The  dispersion  curves  are  an  excellent  fit  to  the  data. 
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Q  inversion  is  linear,  and  takes  little  computer  time. 
Matrices  are  formed  and  decomposed  on  the  first  iteration,  so  the 
best  q  model  ma>  be  quickly  found.  Again  the  command  PLOT  plots  the 
Q  model  while  the  command  PLOT  DATA  plots  the  spectral  ratios  and 
corresponding  data  fit.  The  spectral  ratios  should  look  like  an 
upward  sloping  line  (see  Figure  10).  If  the  data  points  do  not  have 
an  upward  slope,  it  means  the  attenuation  coefficients  decrease  with 
frequency  and  no  Q  structure  will  fit  the  data.  If  this  happens, 
there  is  probably  something  wrong  with  the  inversion,  or  with  the 
data. 

When  the  Q  inversion  is  performed,  the  estimated  moment  and 
of  the  explosion  are  printed  out.  Save  these  numbers,  since  they 
are  not  printed  out  again. 

Since  the  Q  inversion  produces  a  smoothed  s/Q  model,  it  will 
not  reveal  any  true  discontinuities  in  Q  in  the  earth.  In 
particular,  it  will  not  usually  produce  a  very  low  Q  at  the 
surface.  If  high  frequency  data  have  been  used,  it  may  be  necessary 
to  add  a  low  Q  layer  near  the  surface  after  the  inversion  to  reduce 
high  frequency  ringing. 

Once  the  inversion  is  completed  the  final  step  is  to  output  a 
structure  file.  The  structure  file  is  a  formatted  file,  which  is  in 
the  format  of  an  input  file  to  SYNSRF. 

1.2.3  Generation  of  Synthetic  Seismograms  and  Excitation  Functions 

The  third  stage  of  the  path  correction  procedure  is  to  use  the 
program  SYNSRF  to  generate  an  eigenfunction  file,  and  to  make  a 
synthetic  seismogram  as  the  final  test  of  the  procedure. 

SYNSRF  is  divided  into  four  main  sections.  The  first  section 
requires  input  of  a  model  (structure)  file  and  a  choice  of 
frequencies  for  the  calculation  of  phase  velocities.  The  second 
section  calculates  group  velocities,  eigenfunctions  and  related 
parameters.  The  third  section  calculates  a  synthetic  seismogram. 
The  fourth  section  allows  prints  and  plots  of  the  seismogram, 
spectra  and  dispersion  curves. 
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Figure  10.  Q  structure  and  data  fit  obtained  using  a  DF  of  2.5.  The  Q 
model  is  a  smoothed  version  of  the  true  Q  model  (Figure  1) 
while  the  data  is  the  ratio  of  synthetic  (unattenuated)  to 
observed  spectra.  The  moment  is  obtained  simultaneously.  For 
this  example  the  moment  differed  by  2%  from  the  true  moment. 

21 


S-CCBED 


SYNSRF  is  designed  to  be  easy  to  use  and  to  make  it  easy  to 
correct  mistakes.  If  you  run  one  of  the  sections  and  use  an 
incorrect  number,  or  forget  to  output  a  file,  simply  back  up,  make 
the  cnange  and  run  this  section  again.  None  of  the  sections  is  very 
time  consuming  and  all  previous  inputs  are  remembered  (except  for 
output  files,  since  they  are  closed  after  each  section  is  run). 

The  essential  input  in  part  1  of  SYNSRF  is  the  model  file  (the 
output  file  from  INVERT).  The  default  set  of  frequencies  is  100 
frequencies  evenly  spaced  from  .002  to  .2  Hz.  In  most  cases,  this 
will  be  a  good  range  of  frequencies.  For  observations  at  close 
range  (  <  1000  km)  with  an  instrument  that  does  not  filter  out 
higher  frequencies,  it  may  be  necessary  to  use  frequencies  from  .005 
to  .5  Hz  instead.  The  program  is  currently  dimensioned  to  allow  100 
frequencies  and  three  modes.  After  the  model  file  is  read  in  in 
Part  1,  the  program  looks  through  the  shear  velocities  to  estimate 
root  search  parameters  (minimum  and  maximum  phase  velocities  and  a 
step  size).  The  default  values  are  almost  always  adequate, 
especially  for  the  fundamental  mode.  If  the  step  size  chosen  is  too 
large,  a  root  may  be  missed.  This  will  show  up  as  a  glitch  in  the 
dispersion  curve. 

The  command  RUN  causes  dispersion  curves  to  be  calculated  and 
leaves  the  program  in  Part  2.  The  only  inputs  needed  in  Part  2  are 
the  source  depth  and  a  file  name  for  the  eigenfunction  output  file. 
The  default  source  depth  is  1  kilometer.  If  the  structure  being 
used  is  the  path  structure,  then  the  source  depth  is  irrelevant.  It 
is  necessary  to  specify  an  eigenfunction  file  name  here.  This  file 
will  oe  used  in  two  ways.  It  will  be  read  back  in  Part  3  if  a 
separate  source  region  is  used,  and  it  contains  all  of  the 
information  needed  to  construct  the  Green's  function  that  is  the 
final  path  correction. 

The  command  RUN  now  causes  the  eigenfunctions  to  be  generated 
and  leaves  the  program  in  Part  3.  You  can  make  a  plot  of  the 
dispersion  curves  (phase  and  group  velocities)  here  to  compare  with 
the  observations.  Make  this  plot  before  reading  in  a  new  source 
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region  eigenfunction  file,  since  the  velocities  plotted  are  the 
source  region  velocities. 

Several  quantities  must  be  input  here  to  make  a  synthetic 
seismogram.  If  you  are  using  separate  source  and  path  regions,  an 
eigenfunction  file  for  the  source  region  and  the  eigenfunction  file 
for  the  path  must  be  read  in.  An  instrument  response  must  be  read 
in,  preferably  using  the  polynomial  coefficients  that  were  used  in 
TELVEL  (on  a  formatted  rather  than  a  REALIO  file,  see  file  format 
appendix).  The  distance  and  time  window  must  be  specified  here. 

After  the  distance  is  set,  the  time  window  is  estimated  by  the 
program  using  the  group  velocities.  For  the  purpose  of  making  a 
path  correction,  however,  it  may  be  preferable  to  use  the  time 
window  of  the  original  seismogram  to  make  the  comparison  easier. 

The  command  RUN  now  causes  a  seismogram  for  an  explosion  with 
a  of  one  cubic  meter  to  be  generated,  and  leaves  the  program  in 
Part  4.  Then  plots  of  the  seismogram  and  spectrum  may  be  made,  a 

run  summary  may  be  obtained,  and  surface  wave  magnitudes  may  be 

computed. 

The  seismogram  plotted  here  should  look  like  the  original 

seismogram  (see  Figure  11)  after  removal  of  noise  and  multipathing. 
The  amplitude  should  differ  from  the  original  amplitude  by  a  factor 
approximately  equal  to  the  estimated  in  the  inversion  calculation. 

The  last  step  in  the  procedure  is  to  transform  the 
eigenfunction  file  into  a  form  usaole  by  a  moment  tensor  inversion 
program.  A  short  program  EXCITE  turns  the  eigenfunction  file  into  a 
formatted  file  using  the  excitation  function  notation  of  Kanamori 
and  Given  (1981),  as  explained  in  Appendix  4. 


II.  SHAGAN  RIVER-SRO  PATH  CORRECTIONS 


In  this  section  we  describe,  in  brief,  the  processing  of 
explosion  seismograms  recorded  at  seven  SRO  stations.  For  each  path 
there  is  a  description  of  the  analysis  including  an  estimate  of 
confidence  for  each  path  correction,  plus  five  figures.  The  first 
figure  is  a  narrow  band  filter  estimate  of  the  group  velocity  for 
one  seismogram  along  each  path.  The  NBF  estimate  often  shows  up 

problems  such  as  extreme  multipathing,  high  noise  level,  or 

bifurcated  group  velocity  curves  which  may  cause  difficulty  in  the 
data  analysis.  The  second  figure  is  a  plot  of  the  data  fit,  solid 

lines  being  the  calculated  phase  and  group  velocity  curves  from  our 

final  structure  while  letters  U  (group  velocity)  and  C  (phase 
velocity)  are  observed  data  points-  The  third  figure  is  a  plot  of 
shear  velocity  versus  depth  for  the  structure  obtained  by  inversion, 
the  fourth  is  a  plot  of  Q  versus  depth,  while  the  last  figure  is  a 
listing  of  velocities,  densities,  and  Q  obtained  for  each  model. 
The  three  events  processed  are  explosions  at  the  Shagan  River  test 
site  which  occurred  on  2  December  1979  (No.  318),  29  November  1978 
(No.  312),  and  23  Dune  1979  (No.  313).  The  seismograms  are  shown  in 
Figures  12  through  14. 

In  addition  to  obtaining  a  velocity  and  Q  structure  for  each 

path,  we  also  obtain  an  estimate  of  moment  for  each  path.  By 

assuming  a  known  compressional  velocity  and  density  for  the  source 

region  (we  used  a  »  5000  m/sec,  o  »  2100  kg/m  ) ,  the  result  may  be 

2 

expressed  in  terms  of  4':c  using  the  relation  »  4itpa  The 

results  are  listed  in  Table  1. 

The  last  station  is  considerably  less  reliable  than  the  first 
six  for  reasons  that  will  be  explained  with  the  data.  It  is 
interesting  to  note  that  there  seems  to  be  a  correlation  between  the 
estimated  'i'oo  and  the  quality  of  the  data.  This  may  be  an  indication 
of  energy  lost  due  to  scattering  and  non-plane  layered  effects  along 
the  travel  path.  It  does  not  seem  to  result  from  application  of  the 
phase  matched  filter  time  window  since  synthetic  seismograms  made 
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Figure  12.  Seismograms  recorded  at  SRO  stations  for  Shagan  River  i 

explosion  number  318,  December  2,  1979.  Seismograms  were 
processed  to  obtain  path  corrections  for  stations  KAAO, 

SHIO,  ANTO,  CHTO,  KONO,  GRFO  and  MAJO.  Before  processing, 
a  reduced  time  window  was  chosen  to  eliminate  early  arrivals 
and  glitches  and  the  seismograms  were  demeaned  and  detrended. 
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Figure  13.  Seismograms  recorded  at  SRO  stations  for  Shagan  River 

explosion  number  312,  November  29,  1978.  Seismograms  were 
processed  to  obtain  path  corrections  for  stations  CHTO,  MAJO, 
ANTO,  SHIO,  KONO,  and  GRFO. 
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Figure  14.  Seismograms  recorded  at  SRO  stations  for  Shagan  River 
explosion  number  313,  June  23,  1979.  Seismograms  were 
processed  to  obtain  path  corrections  for  stations  CHTO,  MAJO, 
SHIO,  and  KONO.  Processing  of  stations  KAAO  and  ANTO  was 
attempted,  but  excessive  noise  and  multipathing  made  it 
difficult  to  obtain  stable  group  velocity  curves  at  those 
stations. 
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Table  1 
ESTIMATED  4'« 
(MOMENT  ^  6.60  x  10^^) 


Station 

Oi stance 

Event  No.  318 

Event  No.  312 

Event  No. 

KONO 

4370 

15100 

9980 

15000 

SHIO 

2926 

13400 

9760 

13400 

MAOO 

4918 

13000 

5100 

8100 

GRFO 

4700 

11800 

6700 

— 

AN  TO 

3739 

10300 

5440 

— 

CHTO 

3890 

7200 

5000 

6700 

KAAO 

1885  ' 

5700 

— 

29 
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using  the  estimated  values  of  i**  are  not  consistently  lower  in 
amplitude  than  the  observed  seismograms  (see  Table  2).  Of  course, 
the  estimates  of  foo  also  depend  on  the  accuracy  of  the  instrument 
gain.  These  seismograms  were  provided  by  Data  Services  at  the 
Seismic  Data  Analysis  Center  (SOAC)  after  conversion  to  nanometers 
(at  20  seconds). 

In  general,  the  processing  resulted  in  reasonable  velocity 
structures  and  Q  structures  for  most  of  the  stations  and  should 
provide  a  good  estimate  of  the  path  structure  and  attenuation. 
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Table  2 

TIME  DOMAIN  AMPLITUDE  COMPARISON 
station  OBSERVED/SYNTHETIC  AMPLITUDES  (NM) 


Event  No.  318 

Event  No.  312 

Event  No.  313 

KONO 

765/725 

460/479 

736/720 

SHIO 

606/600 

414/437 

576/600 

MAJO 

371/378 

134/148 

248/236 

GRFO 

361/484 

240/275 

— 

ANTO 

192/156 

94/83 

143/  — 

CHTO 

174/249 

132/173 

178/232 

KAAO 

613/305 

— 

282/  — 
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Figures  15.1  through  15.5  are  on  the  following  pages. 


Figure  15.  Path  1:  SHAGAN-KONO 

Distance:  4372  km 
Azimuth:  72.6* 

Instrument:  ASRO-LP 

Events  Processed:  318,  312,  313 

Description:  This  is  an  excellent  station.  A  very 

clean  group  velocity  curve  was  obtained  with  no  evidence 
of  multipathing  or  any  observational  problems. 
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SHAGAN-KOMO  STRUCTURE 
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Figure  15.5  . 
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Figures  16.1  through  16.5  are  on  the  following  pages. 


Figure  16.  Path  2:  SHAGAN-SHIO 

Distance:  2927  km 
Azimuth:  340.8' 

Instrument:  SRO-LP 

Events  Processed:  313,  312,  313 

Description:  This  is  an  excellent  station.  There  is  a 
small  branch  in  the  group  velocity  curve  at  0.065  Hertz, 
but  otherwise  it  is  very  clean.  The  group  velocities, 
seem  stable  and  inversion  results  in  a  reasonable 
structure  with  a  relatively  low  Q. 
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SHAGAN-SHIO  STRUCTURE 
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Figure  16.5. 
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Figures  17.1  through  17.5  are  on  the  following  pages. 


Figure  17.  Path  3:  SHAGAN-MAJO 

Distance:  4911  km 
Azimuth:  307* 

Instrument:  ASRO-LP 

Events  Processed:  313,  312,  313 

Description:  This  is  a  good  station.  There  is  clear 
evidence  of  multipathing,  but  the  separate  arrivals  are 
removed  by  the  phase  matched  filter.  The  group 
velocities  differ  by  0.1  km/sec  at  0.02  Hertz. 
Processing  three  seismograms  and  averaging  helps  to 
obtain  consistent  results. 
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SHAGAW-MAJO  STRUCTURE 
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Figure  17.5. 
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Figures  18.1  through  18.5  are  on  the  following  pages. 


Figure  18.  Path  4:  SHAGAN-GRFO 

Distance:  4699  km 
Azimuth:  62.7* 

Instrument:  SRO-LP 
Events  Processed:  318,  312 

Description:  This  is  a  good,  although  somewhat  noisy, 
station.  The  group  velocity  plot  shows  many  distinct 
arrivals,  but  the  phase  matched  filter  seems  to  separate 
them  well.  The  final  phase  and  group  velocities  are 
quite  consistent  for  the  two  events  processed. 


50 


S-CUBED 


see  I.  OPTIONtll-HClPI  7  >L 
SEC  I.  OPTIONIH'HELPI  7  )1 
caoup  UELS  PICKED) 

IDX  FRED  CRPU  AflPL 


CO 


•V 


j 


1 


. . . . . 

«•  w 


^<\i<vainini<Mmnj<VAiiii«Mni(UfvrunjojruiMmni*vo 


•  •MMrynjr\icjnji^m<*fmn^^«r'vui«i>vr>4/>u*A0t04A«As0i«f*»f<>r>r-r^aoaBOBa9(0aoo 


51 


S-CUBED 


••V:.-.--..-fA!fai;<* 


SHAGAN-GRFO  STRUCTURE 


I 

depth 

THICK 

alpha 

»£TA 

8HO 

OOP 

1 

S.3344M3 

s.a3344a3 

5.377*883 

3.355*883 

2.581*883 

4.278*882 

a 

t.lS7«««4 

S.a334««3 

6.328*883 

3.379*883 

2.596*883 

4.441*882 

3 

1.7Sa«««4 

S.a33*4«3 

6.177*883 

3.467*883 

2.654*883- 

4.381*882 

4 

a.334*9«4 

S.  833 'HMD 

6.398*883 

3.587*883 

2.732*883 

5.346*882 

s 

a.917*«*4 

5.333'»««3 

6.618*883 

3.718*883 

2.811*883 

5.178*882 

6 

3.SM49A4 

S.8334M3 

6.327*883 

3.332*883 

2.391*883 

4.778*882 

7 

4.«84«9«4 

S.3334M3 

7.844*883 

3.354*883 

2.378*883 

4.441*882 

3 

4.679*4«4 

S.949-»««3 

7.267*883 

4.879*883 

3.851*883 

4.279*882 

9 

S.3S84«*4 

6.31S4M3 

7.582*883 

4.211*883 

3.137*883 

4.266*882 

!• 

6.14l«««4 

7.3ea4««3 

7.748*883 

4.349*883 

3.227*883 

4.393*882 

11 

7.83S44S4 

8.346*M3 

7.372*883 

4.475*883 

3.389*883 

4.629*882 

12 

a.a6i«««4 

l.»a44«*4 

8.136*883 

4.567*883 

3.369*883 

4.917*482 

13 

9.a3S**«4 

1.174+444 

8.284*883 

4.685*883 

3.393*883 

5.184*882 

1^ 

l.tSS^MS 

1 . 34S-»««4 

8.172*883 

4.587*883 

3.382*883 

5,373*882 

IS 

t.aia^MS 

1.5411^4 

8.886*883 

4.539*883 

3.358*883 

5.476*882 

IS 

1.3894««5 

1.7S5+444 

8.822*883 

4.583*883 

3.327*883 

5.535*882 

17 

1.SS14MS 

2.3834484 

8.824*883 

4.584*883 

3.323*883 

5.592*882 

18 

t.823«««S 

2.317*884 

8.892*883 

4.542*883 

3.352*483 

5.665*882 

19 

a.«884««5 

3.655*884 

8.183*483 

4.533*483 

3.385*883 

5.739*882 

Figure  18.5 . 
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Figures  19.1  through  19.5  are  on  the  following  pages. 


Figure  19.  Path  5:  SHAGAN-ANTO 

Distance:  3740  km 
Azimuth;  57.2* 

Instrument:  SRO-LP 
Events  Processed:  318,  312 

Description:  This  is  a  difficult,  noisy  station  to 

process.  There  are  numerous  multiple  arrivals  on  the 
group  velocity  plots.  No  group  velocity  curve  could  be 
found  for  Event  313.  Nevertheless,  fairly  consistent 
velocities  were  obtained  for  the  remaining  two  events 
from  0.02  to  0.07  Hertz,  and  when  inverted,  they  produce 
a  reasonable  velocity  and  Q  structure. 
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Figure  20.1  through  20.5  are  on  the  following  pages. 


Figure  20.  Path  6:  SHAGAN-CHTO 

Distance:  3890^km 
Azimuth:  337.1* 

Instrument:  SRO-LP 

Events  Processed:  318,  312,  313 

Description:  This  is  a  rather  noisy  station  with 

multiple  arrivals  and  multipathing  present.  Most 
difficult  is  a  branch  in  the  group  velocity  curve  above 
0.06  Hertz.  This  causes  the  group  velocity  from  the 
phase  matched  filer  to  jump  to  the  other  branches. 
Averaging  velocities  for  three  events  helped  to  obtain  a 
consistent  result. 
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Figures  21.1  through  21.5  are  on  the  following  pages. 


Figure  21.  Path  7:  SHAGAN-KAAO 

Distance:  1884  km 
Azimuth:  22.1* 

Instrument:  ASRO-LP 
Events  Processed:  318 

Description:  This  is  a  very  difficult  station  to 
process  because  of  a  cleanly  bifurcated  group  velocity 
curve  above  0.04  Hertz.  The  split  is  quite  obvious  in 
the  narrow  band  filter,  but  the  phase  matched  filter 
simply  cannot  handle  it,  wanting  to  jump  back  and  forth 
between  the  branches.  To  process  this  seismogram,  we 
selected  the  upper  group  velocity  curve  and  performed  no 
iterations  with  the  phase  matched  filter.  The  final 
result  (synthetic  seismogram  made  using  the  structure) 
is  a  fairly  good  match  to  the  original,  but  we  cannot 
express  much  confidence  in  the  results.  The  same 
problem  will  arise  in  any  further  data  processing  at 
this  station;  so  it  should  probably  be  avoided  if 
possible. 
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III.  CONCLUSIONS  AND  RECOMMENDATIONS 


/Je  have  developed  a  comprehensive  surface  wave  analysis 
package  and  have  used  it  to  obtain  path  corrections  for  paths  from 
the  Shagan  River  explosion  test  sites  to  several  SRO  stations.  This 
analysis  provides: 

1.  An  average  velocity  and  Q  structure  for  each  path. 

2.  A  scalar  moment  estimate  for  each  event  processed. 

3.  A  Green's  function  for  each  path. 

The  Green's  function  may  be  used  co  perform  a  moment  tensor 
inversion  or  to  construct  a  matched  filter  for  events  from  a  given 
area.  In  addition,  the  processing  provides  a  description  of  station 
quality,  by  identifying  paths  with  strong  multipathing,  scattering, 
or  variaole  group  velocities. 

The  procedure  used  to  obtain  a  path  correction  is  to  use 
narrow  band  filtering  and  phase  matched  filtering  to  obtain  phase 
and  group  velocities  from  surface  wave  seismograms,  invert  these  to 
find  the  average  earth  structure  along  the  path,  use  the  structure 
to  generate  a  synthetic  surface  wave  spectrum  and  invert  the 
synthetic/observed  spectral  ratio  for  moment  and  Q  structure.  For 
the  patns  processed  to  date,  we  find  very  consistent  phase  and  group 
velocities  at  most  stations,  which  in  turn  produce  well  constraiited 
shear  velocity  structures.  The  attenuation  measurements  obtained 
using  a  single  station  are  not  as  accurate  as  velocity  measurements; 
however,  we  are  able  to  obtain  reasonable  Q  structures  and 
consistent  moments  in  most  circumstances.  The  method  used  to  obtain 
Q  structures  is  a  new  and  promising  technique.  However,  a 
comprehensive  error  and  trade-off  analysis  has  not  yet  been 
performed.  This  should  be  done  as  part  of  a  continuing  effort  to 
make  further  path  corrections. 

The  interactive  format  of  the  surface  wave  analysis  programs 
allows  data  to  be  processed  very  quickly.  Once  seismograms  and  all 
relevant  information  have  been  obtained  for  a  given  path,  a  path 
correction  can  be  obtained  in  about  two  hours. 
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The  procedure  we  have  developed  works  well  and  is  the  best 
method  available  for  obtaining  path  corrections  from  a  known  test 
site  to  a  set  of  receiving  stations.  There  are  two  aspects  of  the 
problem,  however,  which  have  not  yet  been  addressed  and  which  may  be 
important  to  the  yield  estimation  and  discrimination  problems. 
First,  the  method  can  only  be  applied  to  known  test  areas  so  that 
seismograms  for  several  large  explosions  are  available  for 
processing.  While  the  method  could  be  applied  to  a  single  isolated 
explosion,  it  is  preferable  to  have  more  than  one  event  available  in 

order  to  cneck  the  stability  of  the  results  based  on  internal 

consistency.  Furthermore,  a  moment  tensor  inversion  cannot  be 
performed  until  all  relevant  path  corrections  for  a  given  source 

area  have  been  determined.  Second,  the  method  assumes  that  the 
earth  can  be  adequately  modeled  by  plane-layers  with  constant 

velocities  and  constant  Q.  While  this  approximation  is  usually 
adequate,  the  effects  of  lateral  heterogeneities  and  scattering  are 
often  very  clear  in  the  initial  stage  of  data  processing. 

In  order  to  apply  the  path  correction  procedure  to  new  regions 
of  the  earth,  the  algorithms  could  be  modified  to  allow  the  use  of 
earthquakes  as  well  as  explosions.  A  small  amount  of  code 
development  would  be  required  to  extend  the  method  to  use 
earthquakes.  No  changes  are  required  for  structure  inversion,  but  Q 
inversion  requires  the  construction  of  a  synthetic  spectrum  using  an 
earthquake  (double-couple)  source  instead  of  an  explosion  source. 
Of  course,  application  of  the  method  to  an  earthquake  requires 
knowledge  of  the  earthquake  focal  mechanism  and  depth.  This 
information  can  be  obtained  with  sufficient  accuracy  from  body  wave 
observations.  The  advantage  of  applying  the  method  to  earthquakes 
is  that  earthquakes  occur  almost  everywhere.  If  an  isolated 
explosion  occurs,  then  earthquakes  in  the  area  can  be  used  to 
constrain  the  velocity  and  Q  structures. 

A  crucial  question  in  the  application  of  surface  waves  to 
yield  estimation  is  how  the  surface  waves  are  affected  by  lateral 
structure  variation  and  scattering.  The  use  of  phase  matched 
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filtering  removes  well  separated  multipath  arrivals  and  scattered 
waves  in  order  to  obtain  more  accurate  group  velocities.  It  is  not 
clear  what  effect  these  multiple  arrivals  have  on  the  observed 
amplitudes,  however.  If  there  is  a  large  amount  of  scattered 
energy,  should  it  be  included  in  the  primary  arrival  under  the 
assumption  that  it  was  scattered  from  the  primary  arrival,  or 
excluded  since  we  do  not  know  the  origin  of  the  scattered  waves,  or 
snould  we  avoid  paths  which  show  a  large  amount  of  scattering? 

To  summarize,  we  have  developed  a  comprehensive  surface  wave 
analysis  package  which  allows  path  corrections  to  be  quickly 
determined.  We  recommend  the  following  items  as  topics  for 
continuing  research. 

1.  A  comparison  of  our  observations  of  multiple  arrivals  and 
scattering  with  theoretical  studies  of  these  effects  in 
order  to  estimate  their  influence  on  moment  and 
attenuation  estimates. 

2.  A  comprehensive  error  and  tradeoff  analysis  of  the  Q 
inversion  algorithm. 

3.  Extension  of  the  method  to  allow  the  use  of  earthquakes 
to  obtain  path  corrections  in  isolated  areas. 
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Earth  Structure  File 


Use:  Written  or  read  by  INVERT 

Read  by  SYNSRF  (Part  1  and  Part  2) 

Format:  Read  by  following  statements 

READdEMF,  1)  HED 
1  F0RMAT(A80) 

READdEMF,*)  NLAYER,  NUNITS 
DO  10  I  =  1,  NLAYER 

10  READdEMF,*)  THKd),  AL(I),  BEd),  RHOd),  QM(I),  QB(I) 

Definitions: 

lEMF  ■  File  number 
HED  a  80  Character  header 
NLAYER  a  Number  of  layers 
NUNITS  a  0  for  MKS  units 

a  1  for  Seismological  units  (km,  km/sec,  gm/cm  ) 

THK(I)  a  Thickness  of  layer  I 
AL(I)  a  Compressional  velocity 
BE(I)  a  Shear  velocity 
RHO(I)  a  Density 
QM(I)  a  Shear  quality  factor 
QB(I)  a  Bulk  quality  factor 

Notes:  Free  format  read  is  used  to  make  input  easier.  This 

means  variables  do  not  have  to  be  lined  up  in  columns, 
but  means  that  six  numbers  must  be  input  on  every  line. 
If  either  the  bulk  Q  or  shear  Q  is  set  equal  to  zero,  the 
program  will  set  them  to  a  large  number.  Both  INVERT  and 
SYNSRF  use  MKS  units  internally,  but  seismological  units 
may  be  used  as  input  for  convene! ence. 


Observed  Oispersion  File 


Use:  Written  by  TELVEL 

Also  written  by  AVEDAT,  INTERQ 
Read  by  INVERT,  AVEDAT,  INTERQ 

Format:  Read  by  the  following  statements 

READ(ITEL,1)  HED 

1  F0RMAT(A80) 

READ(ITEL,2)  NF,  GIST,  OLYT,  PHACOR 

2  FORMAT(I5,3E15.5) 

DO  10  I  »  1,  NF 

10  REA0(ITEL,3)  ARL(I),  MODE(I),  F{I),  U{I),  C(I),  DU(I),  DC(I), 
+AMP(I),  OAMP(I) 

3  FORMAT  (A1,I4,7E10.4) 

Definitions: 

HED  —  30  Character  header 

NF  —  Number  of  frequencies 

OIST  --  Distance  from  source  to  receiver  (km) 

DLYT  —  Delay  time  (start  of  original  seismogram) 

PHACOR  —  Initial  source  phase 

ARL  —  Single  character  variable,  R  for  Rayleigh,  L  for  Love 
MODE  —  Mode  number  of  observed  data  point  (usually  1) 

F  —  Frequency  of  observed  data  point 
U  —  Group  velocity 
C  —  Phase  velocity 
DU  —  Standard  deviation  in  U 
DC  —  Standard  deviation  in  C 

AMP  —  Observed,  instrument  corrected  spectral  amplitude 
DAMP  —  Standard  deviation  in  AMP 


Notes: 


The  true  observed  phase  is  given  by 
^  ^  PHACOR 

34 
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Units  are  in  kilometers  for  distance,  km/sec  for  velocities, 
^  transformed  input  units  for  amplitude. 

DU,  DC,  DAMP  are  set  to  zero  by  TELVEL,  and  are  calculated  by 
post  processing  utility  codes  (AVEDAT,  INTERQ) 

For  array  data,  several  files  may  be  combined  into  a  single 
f  file,  and  AMP,  DAMP  replaced  with  attenuation  coefficients  y 

for  multistation  Q  inversion  (INTERQ). 
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Dispersion  File 


Use:  Written  by  SYNSRF  (Part  1) 

Also  written  by  INVERT 
Read  by  SYNSRF  (Part  2) 

Format:  Unformatted  file  read  by  following  statements 

READ(IDISP)  HED 
READ(IDISP)  NRL,  NMODE 
DO  10  I  -  1.  NMODE 

READ(IDISP)  NF,  (F(J,I),  J  =  1,NF),  (C(J,I),  J  =  1,NF) 
10  MF(I)  =.  NF 

Definitions: 

HED  —  80  Character  header 
NRL  —  1  for  Rayleigh 
2  for  Love 

NMODE  —  Number  of  modes 

MF(I)  —  Number  of  frequencies  for  mode  I 
F(J,I)  —  Frequency  J  for  mode  I 
C(J,I)  —  Phase  velocity  0  for  mode  I 

Notes:  F  and  C  are  double  precision  variables  on  this  file 
All  units  are  MKS 
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Eigenfunction  File 

Use:  Written  by  SYNSRF  (Part  2) 

Read  by  SYNSRF  (Part  3)  for  source  or  path 
Read  by  EXCITE 

Format:  Unformatted  file  read  by  following  statements 

READ(IEIG)  HED 

READ(IEIG)  NRL,  NMODE,  DEPTH,  LAYER,  ALPHA,  BETA,  RHO 
DO  10  I  =1  1,  NMODE 

READ(IEIG)  NF,  (F(J,I),  3=1,  NF),  (CU(J,I),  J  =  1,  NF), 
(U(J,I),  J  »  1,  NF),  (AR(J)I),  J  -  1,  NF), 

(EL(J,I),  J  -  1,  NF),  (6A(J,I),  J  =  1,  NF), 
((EIGEN(J,K,I),  J  -  1,  NF),  K  =  1,  NQ) 

10  MF(I)  =  NF 

Definitions: 

HED  --  80  Character  header 
NRL  --  1  for  Rayleigh 
2  for  Love 

NMODE  —  Number  of  modes 

DEPTH  —  Source  depth 

LAYER  —  Source  layer  number 

ALPHA  —  Source  layer  compressional  velocity 

BETA  —  Source  layer  shear  velocity 

RHO  --  Source  layer  density 

MF(I)  —  Number  of  frequencies  for  mode  I 

F(J,I)  —  Frequency  J  for  mode  I 

C(J,I)  ~  Phase  velocity 

U(J,I)  —  Group  velocity 

AR(J,I)  —  Harkrider's  Aj^  (Rayleigh)  or  Aj^  (Love) 

EL(J,I)  --  Ellipticity 

GA(J,I)  —  Gamma,  attenuation  coefficients 

NQ  —  Number  of  eigenfunctions  »  4  for  Rayleigh,  2  for  Love 

EIGEN(J,K,I)  --  Eigenfunction  K  for  mode  I  and  frequency  J 

evaluated  at  source  depth 
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Notes:  All  quantities  are  single  precision 

All  quantities  are  in  MKS  units 
Notation  follows  Harkrider  (1964,  1970) 
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Seismogram  File 


Use:  Written  by  SYNSRF  (Part  3) 

Read  by  RLSEIS  (To  convert  to  REALIO  format) 

Format:  Unformatted  file  read  by  following  statements 

READdSYN)  HEDl,  HED2,  HED3 
READ(ISYN)  NF,  DF,  (SPEC(I),  I  =  1,  NF) 

READ(ISYN)  NT,  DT,  TMIN,  (SEIS{I),  I  =  1,  NT) 

Definitions: 

HEDl,  HED2,  HED3  are  3-80  Character  headers 
NF  —  Number  of  frequencies 
DF  —  Frequency  spacing 

SPEC(I)  —  Complex  spectrum  at  frequency  (I  -1)*0F 
MX  —  Number  of  time  points  in  synthetic  seismogram 
OT  —  Time  spacing 
TMIN  -  Start  time 

S£IS(I)  —  Seismogram,  time  series  at  time  TMIN  +  (I  -1)*0T 

Notes:  All  quantities  are  in  MKS  units 

Spectrum  is  stored  as  real  and  imaginary  parts 
Transform  of  spectrum  will  start  at  T  =  0,  not  T  =  TMIN 
Fourier  transform  is  defined  by: 
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Decomposed  Spectral  File 
Use:  Written  by  SYNSRF  (Part  3) 

Format:  Unformatted  file  read  by  following  statements 

READdSPEC)  HEDl.  HED2.  HED3 

REAO(ISPEC)  NMOOE,  NF,  R,  AMO,  ALPHA,  BETA,  RHO 

READ(ISPEC)  (F(I),  I  .  1,  NF) 

READdSPEC)  (ArNST(J),  PINST(I),  I  -  1,  NF) 

READ(ISPEC)  (ASRC(I),  PSRC(I),  I  -  1,  NF) 

DO  ID  J  -  1,  NMODE 

10  READ(ISPEC)(AMP(J,I),PHASE(J,I),  GAMMA(J,I),  I  -  1,  NF) 
Definitions: 

HEDl,  HED2,  HED3  —  3-80  Character  headers 
NMOOE  —  Number  of  modes 
NF  —  Number  of  frequencies 
R  —  Distance  (meters) 

AMO  —  Source  moment 

ALPHA,  BETA,  RHO  —  Elastic  constants  of  source  region 
F(I)  —  Array  of  frequencies 

AINST(I),  PINST(I)  —  Instrument  amplitude  and  phase  at 
frequency  F(I) 

ASRC(I),  PSRC(I)  —  Source  (RVP)  amplitude  and  phase 
AMP(I,J),  PHASE(I,J),  GAMMA(I,J)  —  Surface  wave  amplitude, 
phase,  and  attenuation  coefficients  for  mode  J,  and  frequency  I 

Notes:  All  quantities  single  precision,  real,  MKS  units; 

frequencies  are  those  used  for  synthetic,  so  there  may  be 
many. This  file  is  intended  for  further  analysis  of 
individual  parts  of  the  spectrum.  The  spectrum  of  the 
complete  complex  seismogram  S(I)  is  given  by; 

S(I)  -  AINSTd)  *  ASRC(I)  *  CEXP(0.,  PINST(I)  +  PSRC(I))* 

5  AMP(I,J)  *  CEXP(-  R  *  GAMMA(I,J),  PHASE(I.J)) 
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Excitation  Function  (Path  Correction)  File 

Use:  Written  by  EXCITE 

Format:  Read  by  following  statements 

REAOdEX,  1)  HED 
READdEX,  2)  ARL,  NF,  R,  AE 
DO  10  I  -  1,  NF 

10  READ  (lEX,  3)  F(I),  C{I),  AR(I),  GAMd),  EPSd),  GREENd) 

20  READ(IEX,  4)  NFl,  DEPTH 
DO  30  I  =.  1,  NFl 

30  READ(IEX,  5)  Fid),  ARld).  RNd),  RPd),  RQd),  RSd) 

Return  to  20  for  more  depths  until  NFl  =  0 

1  FORMAT(AaO) 

2  F0RMAT(A5,  15,  2E10.4) 

3  F0RMAT(7E10.4) 

4  F0RMAT(I5,  E10,4) 

5  F0RMAT(7E10.4) 

Definitions:  HED  —  80  Character  header 

ARL  —  R  for  Rayleigh 
NF  —  Number  of  frequencies 
R  —  Distance  in  meters 
AE  —  Earth's  radius  in  meters 
F(I)  —  Frequency 

C(I)  —  Phase  velocity  at  frequency  F(I) 

AR(I)  —  Harkrider's  Aj^ 

GAM(I)  —  Gamma,  attenuation  coefficients 
EPS(I)  —  Ellipticity 
GREEN(I)  —  Complex  Green's  function 
DEPTH  —  Source  depth 

NFl  —  Number  of  frequencies  for  this  depth 

F1(I)  —  Frequencies  for  this  depth 

C1(I)  —  Phase  velocity 

ARl(I)  —  Harkrider's  Aj^  at  frequencies  FI 

RN(I),  RP(I),  RQ(I),  RS(I)  —  Excitation  functions  as 

defined  by  Kanamori  and  Given 
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Notes: 


This  file  is  the  interface  between  the  Surface  Wave 
Analysis  Package  and  Sierra  Geophysics  Moment  Tensor 
Inversion  Program. 

This  file  is  constructed  from  any  number  of  Eigenfunction 
files  generated  by  SYNSRF,  so  it  may  or  may  not  have  the 
same  set  of  frequencies  for  each  depth. 

The  Green's  function  GREEN(I)  is  defined  by: 

It  -ittir 

GREEN  -  e^^  e  ^  e'^*" 
sin(r/ag) 

The  radius  of  the  earth  a^  is  written  on  the  file 

e 

because  it  is  used  in  constructing  the  excitation 
functions.  Kanamori  and  Given's  excitation  functions  are 
called  Nj^,  P 1^,  Qq,  S|^. 
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Use:  Generated  by  S-CUBED  explosion  codes 
Read  by  SYNSRF 

Format:  Read  by  following  statements 

READ(IRVP,1)  HED 
REA0(IRVP,1)  HED2 
READ{IRVP,2)  YIELD 
REA0(IRVP,1)  JUNK 
READ(IRVP,3)  REL,  VP,  VS 
10  REA0(IRVP,4)  F{I),  AMRVP(I),  PHRVP(I),  LAST 
AM.^VP(I)  -  AMRVP(I)  *  l.E-6 
IF(LAST.  NE.l)  Go  to  10 

1  F0RMAT(A80) 

2  FORMAT (lOX,  E10.4) 

3  F0RMAT(3F10.2) 

4  FORMAT(E15.7,  30X,  2E15.7,  15) 

Definition: 

HEDl,  HED2  —  80  Character  headers 

YIELD  —  Original  explosion  yield  for  calculation 

JUNK  —  Nothing,  blank  line 

REL  —  Elastic  radius 

VP,  VS  —  Elastic  constants  of  source 

F(I)  —  Frequency 

AMRVP(I)  —  Amplitude  of  RVP,  scaled  to  meters  cubed 
PHRVP(I)  —  Phase  of  RVP 

Notes: 

Check  with  Norton  Rimer  for  accuracy  of  YIELD,  VP,  VS. 
To  scale  from  yield  YOLD  to  YNEW 

AMRVP  -  YNEW/ YOLD 
F  -  (YOLD/YNEW)  **  (1/3) 

REL  -  (YNEW/YOLD)  **  (1/3) 


Instrument  Files 


Use:  Read  by  SYNSRF 

Types:  There  are  three  types  —  polynomial  coefficients 

(preferred),  and  two  formats  for  tabular  responses, 
corresponding  historically  to  S-CUBED  long  period  and 
short  period  instrument  responses. 

Format  of  Polynomial  File: 

The  following  format  is  used  for  polynomial  files: 

R£AD(INST,1)  HED 
READ{INST,*)  INUM 
IR  *  INUM  +  1 

REA0(INST,*)  (POLYN(I),  I  -  1,  IR) 

READ(INST,*)  IDEN 
IR  -  IDEN  +  1 

REA0(INST,*)  (POLYO(I),  I  «  1,  IR) 

Definitions:  A  polynomial  response  is  written  in  the  form 


P(s)  -  a.s' 

"T^ - : 

i-'o  Si^' 


where  s  *  iw 

N  a  INUM  the  degree  of  the  numerator  polynomial,  with 
coefficients  »  POLYN(I)  starting  from  degree  zero 
M  a  IDEN  the  degree  of  the  denominator,  with  coefficients 
Bj  a  POLYD(I),  starting  from  degree  zero 

Note:  A  program  (ZERPOL)  is  available  to  transform  poles  and  zeros 
into  polynomial  coefficients.  Another  (RLINST)  puts  a 
polynomial  file  into  REALIO  format  for  use  with  TELVEL  and 
other  REALIO  routines. 
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Format  of  S3LP  file: 


REA0(INST,1)  NF  . 

00  10  I  -1,  NF 

READ(INST,2)  PERIOD,  AMP(I),  PHM2PI 
FR(I)  =  1. /PERIOD 
10  PHS(I)  -  -  2.  *  PI  *  PHM2PI 

1  F0RWT(I3) 

2  F0RMAT(F8.0.  2F12.0) 

Definitions: 

NF  »  Number  of  frequencies 
AMP  »  Instrument  amplitude 
PHS  -  Instrument  phase 
FR  «  Frequencies 

Note:  Phase  stored  as  minus  phase/ 2ir 

Format  of  S3SP  file: 

READ(INST,1)  HED 
'read (INST, 2)  NF 
DO  10  I  -  1,  NF 

10  REA0(INST,3)  AMP{I).  PHS(I),  FR(I) 

1  F0RMAT(A80) 

2  FORMAT (I 10) 

3  F0RMAT(3E15.5) 

Notes:  The  S3LP  and  S3SP  notations  are  historical  only.  An 

instrument  response  may  be  used  with  either  format, 
regardless  of  whether  it  is  a  long  period  or  short  period 
instrument. 

Warning:  Some  instrument  files  have  a  factor  of  minus  one  built 

in,  apparently  to  change  vertical  down  to  vertical  up  in 
old  surface  wave  codes.  Check  before  using. 


REAL  10  Files 


REAL  10  files  are  designed  for  internal  use  by  subroutines  in 
the  time  series  analysis  package  (Berger,  ^  ,  1979)  which  are 

used  by  TELVEL.  The  exact  format  of  the  REALIO  files  is  machine 
dependent,  so  it  is  not  given  here,  however  files  may  be  read  or 
written  by  a  call  to  suoroutine  REALIO.  The  following  commands 
cause  a  file  to  be  opened  and  the  array  SEIS  containing  NPTS  points 
read  in.  The  underscore  is  used  as  a  delimiter  here. 

CHARACTER*80  FILENM 

DIMENSION  SEIS  (NPTS) 

CALL  GETSTR  ('ENTER  SEISMOGRAM  FILENAME  ,  FILENM) 

CALL  GETLU  (FILEfW,  LU,-1) 

CALL  REALI0(SEIS,1,NPTS,LU,0) 

The  call  to  GETSTR  reads  the  file  name  from  the  keyboard. 

Tne  call  to  GETLU  returns  a  logical  unit  number  and  opens  the 
file.  The  -1  used  as  the  third  parameter  indicates  an  existing 
file.  A  number  greater  than  zero  opens  a  new  file. 

The  call  to  REALIO  reads  the  file  into  the  array  SEIS.  The 
use  of  zero  as  the  fifth  parameter  causes  the  file  to  be  read. 
Similarly  if  this  parameter  is  greater  than  zero,  the  array  is 
written  to  tne  file. 

For  more  information  on  the  use  of  these  routines,  see  the 
report  by  Berger  e^  al_. ,  (1979).  The  program  RLSEIS2  is  easily 
modified  to  put  any  file  into  REALIO  format. 
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Several  short  programs  are  available  for  use  with  the  surface 
wave  analysis  package.  The  most  important  are  listed  here. 

EXCITE  —  Converts  eigenfunction  files  to  an  excitation 
function  file.  Requires  only  file  names  as  input. 

RLSEIS  —  Converts  seismogram  file  from  SYNSRF  to  REALIO 
format . 

RLSEIS2  —  Reads  in  a  seismogram,  plots,  truncates,  demeans, 
detrends,  and  tapers,  and  writes  file  in  REAlIO  format. 
Easily  modified  to  read  any  format.  Requires  file  names  and 
number  of  points  as  input. 

RLINST  —  Converts  polynomial  instrument  file  to  REALIO  format. 

AVEDAT  —  Reads  TELVEL  files,  finds  average  phase  and  group 
velocities  and  standard  deviations,  and  relative  moments. 
Outputs  a  new  TELVEL  file  for  path. 

INTERQ  —  Reads  TELVEL  files  for  an  array  of  stations,  finds 
attenuation  coefficients  and  average  phase  and  group 
velocities  and  standard  deviations.  Outputs  a  new  TELVEL  file 
with  amplitudes  replaced  by  attenuation  coefficients. 

ZERPOL  —  Turns  list  of  poles  and  zeroes  for  instrument 
response  into  polynomial  coefficients. 
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APPENDIX  3 

surface  wave  notation 
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The  surface  wave  code  SYNSRF  uses  the  notation  of  Harkrider 
(1964,  1970)  for  surface  waves  with  some  minor  modifications 

displacement  is  given  by: 

X(9,h)A-  >.i(l‘^2m)T  -ikr^-Yr 

u  »  '  '  '  R  S((i))e  e  e  e 

\/2ir<i)CagSin(r/ag) 

where  x  are  functions  defined  by  Harkrider  (1970),  S(u)  is  a  source 
function  and  m  -  1  for  double  couple  and  0  for  point  force  or 
explosion. 

e  -  -1  for  vertical  component  Rayleigh  wave  (positive  up) 
e  ■  ie,  where  e  »  ellipticity  for  radial  component 
e  a  i  for  Love  wave 

For  an  explosion: 

S  a  -3iru  fao 
X  a  (E^  -  E^/Eu) 

For  an  impulsive  point  force: 

S  a  ilgC  (Ig  a  Impulse) 

X  a  -^2  (downward) 

X  a  iEj^  cose  (radial ) 

For  a  double  couple: 

S  a  iMg  (*^0  “ 

X-  C(2fl^/a^  -  3/2)E^  -  E3/pa^]  ^  COS  2e  45*  dip  slip 

X-  -Ej^  sin  2e  strike  slip 

Here  e  is  clockwise  from  north  a  -  azimuth  which  is  used  in 
interactive  input.  Note  that  near  the  free  surface  at  long  periods 
E^  <  0,  E2  >  0,  and  e  <  0,  while  Ej/u  is  small.  Thus  the 
initial  phase  may  be  found  from  these  expressions.  For  an 
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explosion 


<^0  *  -  T 

For  a  downward  point  force 


For  a  point  force  where  velocity  is  computed,  multiply  by  iw  so 

A  3 
^0  ’  ■  T  ’• 

We  have  used  the  following  shorthand  notation  for  Harkrider's 
eigenfunctions: 

-  U*(h)/WQ 

E2  »  w(h)/wQ 
£3  -  a*(h)/(wQ/c) 

£4  =■  t(h)/(wg/c) 
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APPENDIX  4 

CONVERSION  OF  HARKRIDER  EIGENFUNCTIONS  TO 
KANAMORI  AND  GIVEN  EXCITATION  FUNCTIONS 
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We  want  to  convert  the  eigenfunctions  defined  by  Harkrider 
(1964,  1970)  into  the  excitation  functions  defined  by  Kanamori  and 
Given  (1981).  The  easiest  way  to  do  this  is  by  comparing  their 
final  expressions  for  farfield  Rayleigh  waves.  In  the  notation  of 
Harkrider,  the  vertical  component  of  the  Rayleigh  wave  (positive  up) 
is  given  by: 


At,  ■— 

u  ,  -  3  ,  ..  x(a.h)  e~^^  e  (1) 

y/  ZirwcagSina 

The  function  X  depends  on  source  type  and  source  depth.  For  an 
explosion  with  a  step  function  time  history: 

X  »  Siru'i'^e  ~ 

And  for  a  double  couple  with  a  step  function  time  history: 

X  ■  e  (“^Q  iCd^sine  +  d^coss]  sin2©  d^cosZa)  (3) 

where  the  functions  d^  are  defined  by  Harkrider  (1970,  page 
1938).  Here  »  is  measured  west  from  north. 

In  the  notation  of  Kanamori  and  Given,  using  a  moment  tensor 
representation  for  the  source,  again  assuming  a  step  function  time 
history,  the  farfield  vertical  displacement  is  given  by: 


u 


F(®,h) 


(4) 


107 


S-C^BED 


where 


ns.h)  -  -  Sin2»  -  ^  COS29]  (5) 

"  7  [Sr  '  "rI 

*  5  -  Sa](»„  •  ",y) 

The  equivalence  between  the  excitation  function  Pj^,  Sj^, 

N  and  the  Harkrider  eigenfunctions  may  be  found  by  comparing 
equations  (1)  and  (4)  for  special  cases.  For  an  explosion, 
^xx  "  \y  “  ^zz  “'^o’  eRuation  (5)  reduces  to 

F(s.h)  -  M^Nj^  (6) 

The  moment  of  an  explosion  is  defined  by 

a  4iroa^'i'„  -  4ir(A  ■*■  2u}'¥„  (7) 

SO  comparing  equations  (1)  and  (4) 

(8) 

Second,  consider  a  strike  slip  double  couple.  The  only 
nonvanishing  coefficient  in  Equation  3  is  d^  »  -Ej^.  Using 
^xy  “  "^o‘  Equation  (5)  reduces  to  F  -  +  M^P,^  sin2e. 
Comparing  equations  again,  we  get 

(9) 
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r~ 


f 


9 


9 


« 


Next,  consider  a  horizontal  dip  slip  double  couple 

have  d2  ■  £^/w  and  F  *  iM^Qj^  cos©. 
So 


u 


zz 


Fi nal ly, 

M  . 

O’  yy 


a  45* 
-  and 


dip 


slip 


double 


(10) 


couple  has 


1  rax 

?L“ 


2u 


T7‘'3 


Me  also  have 


[- 


1 

I 


cos2©  + 


1 

I 


J  0 


So 


2u  p  4.  . .2  c  "! 

^  -y/ZituCag 

L  T  +  2u  H  X  +  2u  ^3  J 

(11) 


Note  that  these  excitation  frunctions  differ  by  factors  of  -1, 
i(tt,  and  u  from  the  excitation  functions  defined  by  Ben-Menahem, 
Rosenman  and  Harkrider  (1970).  This  is  a  result  of  the  delta 
function  time  history,  and  source  orientation  defined  in  terms  of 
fault  plane  parameters  rather  than  moment  tensor. 


t 


X 


I 
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Bache,  Rodi  and  Harkrider  (1978),  describe  a  method  for 
synthesizing  Rayleigh  wave  seismograms  when  the  source  region 
structure  and  propagation  path  structure  differ.  In  the  following 
expressions  the  Index  1  refers  to  the  source  region,  while  2  refers 
to  the  propagation  path.  The  vertical  component  of  the  Rayleigh 
wave  is  given  by: 


yj  ZirwC^agSina 


X(©,h) 


(1) 


The  radial  component  is  given  by 


u  *  -ic2U 

where  62  is  the  ellipticity  in  region  2. 

The  excitation  functions  described  in  the  last  appendix  can  be 
modified  for  separate  source  and  travel  path  structures,  by 
computing  them  according  to  Equations  (8)  through  (11)  of  the  last 
appendix  using  the  source  region  structure,  and  then  multiplying  by 


Explicitly,  the  vertical  component  of  the  Rayleigh  wave  is 
given  by: 


U  “  F(Nj^,  P|^,  Qq,  Sj^)  X  X  ^ 


TT  -luir 

iT  W 

e  X  e  X  e 


VSin 


where  a  ■  with  a^  the  earth's  radius,  and  F  the  function 

defined  in  Equation  5  of  the  last  section.  To  generate  surface 
waves  from  a  source  described  by  a  moment  tensor,  the  following 
quantities  are  needed: 
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Source  region: 

Nr,  Pr,  Qr,  Sr,  C^, 

Path  region: 

^2’  ^R2’  ''2’  ®2 

In  addition,  the  set  of  frequencies  for  each  of  the  above  functions 
must  be  known,  and  the  distance  r  and  earth's  radius  must  be 
specified.  For  the  earth's  radius,  we  use  a^  »  6.371  x  10^  m. 
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APPENDIX  6 

INVERSION  FOR  Q  AND  MOMENT 


♦ 
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The  surface  wave  displacement  spectral  amplitude  is  given  by: 

u  -  S^(w)  (1) 

2 

where  is  a  constant  -  4itpa  '■?„  for  an  explosion,  S^(u)) 

depends  only  on  the  source  structure  and  $2(11))  and  yI'^*)  depend 
only  on  the  path  structure.  Explicitly: 


(2) 


S2(<-) 


(3) 


where  and  frequency  dependent  structure  excitation 


functions  for  source  and  path  regions,  Cj^ 

rce  and  path  regions, 

,6 


and  C2  are  phase 
velocities  for  source  and  path  regions,  a.  is  the  earth's 

g  “ 

radius  -  6,371  x  10  M,  r  is  the  observation  distance  and 
a  »  r/a  .  and  E^^  are  Harkrider  eigenfunctions  and  0,  a,  s 

are  source  layer  density  and  velocities. 


Note  that  Sj^(id)  is  known  for  an  assumed  source  region 
structure,  and  S2(<>t)  is  computed  for  the  inverted  structure  within 
the  inverse  code.  In  Equation  (1)  if  u  is  the  observed  spectrum, 
then  only  and  y(‘»)  are  unknowns.  Taking  the  logarithm  of 
Equation  1, 


£n(u((o))  -  fin(M^)  +  ln{S^{,a)  )  -  »' 

or: 


t(w) 


£n(M^) 


(Si(«)S2(u.)) 

uTIH 


(4) 


(5) 
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But  we  also  know  the  relationship  between  v  and  Q  (e.g.  Mitchell, 
1975) 


'(to) 


U) 

I? 


ac 

3B 


8 

i  "3 


t 


3C 


*£  1 

5“ 


(6) 


!  *e  I 


where  I  represents  a  sum  over  layers  in  the  earth  model, 


is  not  independent  of  Q 


then 


8. 


If  bulk  losses  are  small, 


So  we  can  write 

y(<u)  *  G((j),h)  (8) 

3 

where  G  is  the  operator  defined  in  Equation  (6),  or  in  matrix  form: 


G  ( (j ,  h ) , 


11 


3 

r 

»  2n 

_£nMo_ 

S^{iu)S2(u))  j 

- ^ - J 


(9) 


These  equations  are  solved  by  INVERT.  In  this  way,  the  Q  structure 
and  Moment  may  be  obtained  simultaneously.  For  array  data.  Equation 
(8)  is  solved  by  INVERT. 


APPENDIX  7 
PROGRAM  TELVEL 
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This  appendix  describes  in  detail  the  logic  of  TELVEL  and  the 
nteractive  input/output  quantities  needed  in  its  execution. 

TELVEL 

A  BRIEF  DESCRIPTION  OF  THE  MAIN  SUBROUTINES 

DEBLIP  -  Removes  blips  from  input  seismograms. 

DETRND  -  Removes  the  mean  and/or  linear  trend  from  a  given 
series. 

FFTR  -  Fast  Fourier  transform. 

FIL  -  Narrow  band  filter  Gaussian  filter  normalized  so  that 
the  integral  over  frequency  from  minus  infinity  to 
infinity  is  equal  to  one.  Explicitly 

FIL(F,  FC,Q)  -VT  (- 

where 
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FUNRAP 


GETTS 


GPLOT 


GPL0T2 


GRPOLY 


Phase  unwraps  the  phase  function  in  the  calling 
arguments.  The  user  specifies  a  threshold,  defining  a 
2it  discontinuity.  Generally  the  threshold  should  be 
between  2  and  3. 

Retrieves  a  segment  of  data  from  the  input  device. 
This  version  assumes  a  REALIO  file.  Segments  of 
length  2**N  returned. 

Plots  group  velocity  curves  from  the  TPK  and  APK 
vectors  generated  by  N8F.  GPLOT  also  plots  the  “best" 
group  velocity  curve  as  determined  by  GRPV.  The 
"best"  curve  will  be  plotted  as  a  dotted  line,  the 
other  peaks  as  integers  representing  the  peak 
amplitude  as  a  fraction  of  the  largest  peak  in  the 
plot  (in  tenths) . 

Plots  old  and  new  group  velocity  curves  to  enable  the 
user  to  make  the  decision  as  to  whether  another 
iteration  is  necessary. 

Estimates  the  group  delay  in  a  system  function  by 
approximating  the  phase  derivative  as 


The  phase  derivative  is  estimated  by  evaluating  the 
phase  of  the  system  function  at  points  on  either  side 
of  the  frequency  of  interest.  If  the  first 
differences  are  equal  within  a  certain  tolerance,  the 
phase  derivative  is  computed  by  taking 

U 
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otherwise,  the  step  size  is  decreased  and  the  process 
repeated . 

GRPV  -  Computes  a  group  travel  time-frequency  file  for  a 
section  of  a  seismogram  and  plots  the  di sponsion 

curve.  With  GVFIND,  GRPV  will  also  find  the  ''best'* 

dispersion  curve. 

GVFIND  -  Picks  the  dominant  group  velocity  curve  from  a  group 
travel  time-frequency  file.  To  be  accepted,  a  peak 
must  be  within  a  certain  time  window  near  the  arrival 
time  of  a  peak  in  an  adjacent  frequency  band.  In 

addition,  the  peak  must  be  above  a  specified  amplitude 
threshold.  If  several  peaks  from  one  filter  pass 

these  criteria,  then  an  error  “radius"  is  computed  for 

each  peak  based  on  its  arrival  time  and  amplitude. 

The  peak  with  the  smallest  radius  will  be  accepted. 

GVITER  -  Handles  the  iterative  portion  of  the  velocity 

determination.  GVITER  assumes  the  user  will  always 

want  to  compare  results  from  the  previous  iteration 
with  the  results  from  the  current  one  before  making  a 
decision  as  to  whether  another  is  necessary.  Thus, 
the  first  step  in  GVITER  is  always  to  compute  new 

group  velocity  curves  from  the  phase  spectrum 

resulting  from  the  previous  iteration.  Steps  are  as 
follows: 

1)  in  a  loop,  for  each  section 

a)  get  from  the  ^  file 

b)  compute  d^/du 

c)  compute  equivalent  group  velocity  and  put 
answer  in  temporary  arrays 
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2)  ask  the  user  which  sections  he  wants  to  look 
at  and  plot  both  old  and  new  group  velocity 
curves 

3)  ask  the  user  if  he  wants  a  new  iteration: 

a)  yes  -  transfer  group  velocities  from 
temporary  arrays  into  GV  and  FV,  then 
exit.  The  main  program  will  continue 
with  PHAV 

b)  no  -  ask  the  user  if  he  wants  the  new  or 
old  GV's,  then  exit  and  the  main  program 
will  go  to  POUT 

4)  upon  return,  ITER  »  1  if  a  new  iteration  was 
requested.  0  a  no. 

IiNTERP  -  Finds  the  estimated  exact  location  and  amplitude  of 
the  peak  or  trough  by  fitting  a  parabola  to  three 
adjacent  trace  points. 

LOCPCK  -  Locates  peaks  in  the  filter  envelope  function.  LOCPCK 
attempts  to  do  the  location  using  the  three  point  rule: 

Xl(I-l)  <  X1(I) 

X1(I+1)  >  X1(I) 

All  resulting  peak  times  and  amplitudes  are  saved. 
The  total  number  of  peaks  saved  is  returned  in  NTOT. 

MPLOT  -  Plots  the  matched  filters  generated  in  PHAV.  The 
Fourier  transform  of  the  matched  filter  is  read  off 
the  REALIO  file,  transformed  to  the  time  domain  and 
plotted.  Only  the  first  LSEC  points  are  plotted  so 
the  filters  can  be  compared  directly  with  the 
seismograms. 


r 


t 


f 


1 


{. 


NBF  -  Performs  all  narrow  band  filter  operations 


1)  FFT  entire  time  series  in  place 

2)  loop  over  filters  and 

a)  produce  filter  output  spectrum 

b)  produce  quadrature  spectrum 

c)  transform  both  to  time  domain 

d)  combine  these  to  produce  envelope  function 

e)  calculate  the  instantaneous  phase 
(optional ) 

f)  locate  the  peaks  in  the  envelope  function 

g)  sort  peaks  in  descending  order 


PHAOE  -  Computes  phase  delay  curve  given  the  arrival  time  of 
peaks  in  the  surface  wave  train.  The  method  used  is 
to  compute  approximations  to  the  group  arrival  time 


TG 


d  phase 
du 


using  the  arrival  time  and  the  period  of  the  peaks. 
The  derivative  is  then  integrated  to  give  the  PHAS 
array.  In  addition,  frequency  -  group  velocity  pairs 
may  be  read  off  an  NBF  -  type  dispersion  curve. 


PHAV  -  Loops  through  the  sections  of  the  seismogram  file  and 
computes  a  phase  delay  curve  for  each  section.  The 
method  is  to  integrate  the  phase  derivative  as  given 
by  the  group  velocity  dispersion  curve  defined  by  the 
vectors  GV  and  FV,  PHAV  writes  out  the  cross 
correlations  to  a  second  and  the  phase  delay  function 
to  a  third.  The  first  is  given  in  its  frequency 

domain  counterparts  since  the  cross  correlating  is 
done  in  the  frequency  domain.  Each  section  of  these 
two  files  contains  2**L0G2N  complex  numbers 

corresponding  to  the  FFT  of  the  associated  time  series. 
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POUT 

PROC 

PVPLQT 

SPLQT 

SYSFUN 


TAPER 


To  perform  the  cross  correlation  in  the  frequency 
domain,  the  time  series  are  extended  with  zeros  to 
length  2**(L0G2N+1) ,  where  L0G2N  is  l+IFIX(log2 
LSEC).  The  extra  zeros  are  necessary  to  avoid 
problems  with  circular  correlation.  Consequently,  the 
output  time  series  will  be  at  least  twice  as  long  as 
the  length  of  each  section. 

Prints  out  a  summary  of  the  run  and  the  phase  and 
group  velocity  curves  as  determined. 

POP  element  with  a  number  of  useful  comments. 

Plots  the  phase  velocities  as  computed  from  the  phase 
functions  in  the  file  referenced  by  LUPHI.  The  user 
has  the  option  of  specifying  the  section  for  which  to 
compute  the  phase  velocity  and  the  number  of  2ir's  to 
add  or  subtract. 

Plots  multiple  seismograms. 

Evaluates  a  system  function  given  in  terms  of  a 
complex  polynomial.  SYSFUN  returns  the  amplitude  and 
phase  of  the  system  function  at  a  specified  frequency. 

If  POLY(l)  <  0,  the  rest  of  the  polynomial  will 
be  ignored  and  the  return  values  will  be 

AMP  -  1. 

PHA  «  0. 

Applies  a  cosine  taper  to  the  ends  of  a  given  series 
and  returns  a  tapered  series. 


126 


S-COBED 


p 


9 


t 


t 


TSEDIT  -  Performs  all  the  time  series  editing  which  may  include: 

1)  deblip  (remove  glitches  from  data) 

2)  demean  (remove  the  mean) 

3)  detrend  (remove  a  linear  trend) 

4)  taper  (apply  a  cosine-bell  taper  to  both  ends) 

5)  insert  leading  zeros  and  update  start  time 

6)  find  the  nearest  power  of  2  for  the  FFT  and  pad 
the  time  series 

7)  some  processed  time  series  on  disc 
VELS  -  The  main  program 

CPLOT  -  Plot  the  cross-correlation  functions  generated  in  PHAV. 
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TELVEL 


SUBROUTINE 


FFTC 

FFTR 

GPLOT 

GRPOLY 

GRPV 

GVFINO 

GVITER 

LOCPCK 

MPLOT 

N6F 

PHAOE 

PHAV 

POUT 

TSEDIT 

VELS 

XPLOT 


CALLS 


COOL 

COOL 

GRPOLY 

SYSFUN 

GETTS,  GPLOT,  GPLT(GPLOT),  GVFINO,  N6F,  TSEDIT 

GRPOLY 

GPLOT2 

INTERP 

FFTR,  PHAOE 

FFTC,  FFTR,  FIL  (function),  LOCPCK 
GO ELAY 

FFTR,  FUNRAP,  MPLOT,  PHAOE,  PPLOT(PVPLOT) ,  PVPLOT, 

SYSFUN,  XPLOT 

SYSFUN 

(OEBLIP),  OETRND,  FFTPRM,  TAPER 

GRPV,  GVITER,  PHAV,  PHITER  (PHAV),  POUT,  SPLOT 

FFTR 


(A)  A  IS  COMMENTED  OUT 
A(B)  A  IS  AN  ENTRY  POINT  OF  B 
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TELVEL 


SUBROUTINE  IS  CALLED  BY 


COOL 

FFTC,  FFTR 

DEBLIJ- 

TSEDIT 

DETRND 

TSEDIT 

FFTC 

NBF 

FFTPRM 

TSEDIT 

FFTR 

MPLOT,  NBF, 

PHAV,  XPLOT 

FIL 

NBF 

F UNRAP 

PHAV 

GDELAY 

PHADE 

GDELAY (GRPOLY) 

GETTS 

GRPV 

GPLOT 

GRPV 

GPL0T2 

GVITER 

GPLT 

GRPV 

GPLT(GPLOT) 

GRPOLY 

GPLOT,  GVFIND 

GRPV 

VELS 

GVFIND 

GRPV 

GVITER 

VELS 

INTERP 

LOCPCK 

LOCPCK 

NBF 

MPLOT 

PHAV 

NBr 

GRPV 

PHAOE 

MPLOT,  PHAV 

PHAV 

VELS 

PH ITER 

VELS 

PHITER(PHAV) 

POUT 

VELS 

PPLOT 

PHAV 

PPLOT(PVPLOT) 

PVPLOT 

PHAV 
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SPLOT 

SYSFUN 

taper 

TSEDIT 

VELS 

XPLOT 


VELS 

GRPOLY,  PHAV,  POUT 

TSEDIT 

GRPV 

PHAV 


A(B)  A  IS  AN  ENTRY  POINT  B 
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TELVEL 


I 


SUBROUTINE  HAS  THE  ENTRY  POINTS 


GPLOT 

GPLT 

GRPOLY 

GO ELAY 

GRPV 

OON8F 

GV ITER 

TNSFER  , 

TRIGSM 

PHAV 

PHITER, 

WINDOW 

POUT 

TRIAVG 

PVPLOT 

CONVRT, 

PPLOT 
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The  following  pages  are  designed  to  enable  the  user  to 
execute  TELVEL  and  also  be  able  to  locate  the  section  of  the  code  in 
which  control  is  currently  assigned.  The  format  of  each  prompt  is 
as  follows: 


^number  of  this  prompt 

f 

#  SUBROUTINE  in  which  prompt  is  located 
"PROMPT  MESSAGE" 

options,  comments,  etc. 
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SEISMOGRAM  FILE  (Q  TO  QUIT) 


Q  -  Stop  execution 

name  -  Name  of  REALIO  seismogram.  Must  have 
constant  at  between  points.  At  end  of 
program,  control  is  restored  to  this  point. 


VELS 

"N  SECTS,  LEN  SECT,  DELTA  -  T,  DELAY  TIME  (+  -  LATE)  “ 


8  - 
N  SECT  - 
LEN  SECT  - 
DELTA  T  - 
DELAY  TIME  - 


Go  back  to  1 

Number  of  sections  in  seismogram. 

Length  of  this  section. 

Time  interval  between  data  points. 

Number  of  seconds  between  the  event  time  and 
the  start  time  of  the  seismogram.  A  plus 
value  indicates  that  the  seismogram  starts 
after  the  event  time  and  a  negative  before. 


VELS 

"DISTANCES  (KM),  OR  NEAREST  OFFSET,  -  INCREMENT  " 

B  -  If  the  first  character  entered  is  "B",  go 
back  to  2. 

distances  -  There  must  be  at  least  two  distances 
specified,  even  though  only  one  section  is 
used.  These  distances  can  be 
input  in  two  ways: 
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a)  -  the  first  value  is  a  distance  in  kai  and 
the  second  value  is  a  delta  distance,  also 
in  km.  This  is  indicated  by  a  negative 
sign  preceeding  the  delta.  The  sign  of 
the  distance  (first  value)  indicates 
whether  the  delta  is  added  to  (positive) 
or  subtracted  from  (negative)  the  starting 
value.  For  this  case  the  algorithm  is: 

DIST(l)  «  ABS(first  value) 

0IST(2)  -  DIST(l)  *  increment 

DIST(3)  -  DIST(l)  ±  2*  increment 

b)  The  actual  distances  in  km  are  entered. 
There  must  be  NSEC  values  and  if  a  “B"  is 
encountered  after  specification  of  the 

second  distance,  control  is  transferred 
back  to  3  and  the  distance  input  starts 
over. 

4  SPLOT 

"SECTIONS  TO  PLOT  -  S  =  ENDLIST,  A  .  ALL,  G  -  GO" 

SPLOT  plots  the  seismogram  sections  requested. 

G  -  Go  back  without  plotting 
B  -  Don't  plot.  Go  back  to  3. 

A  -  Plot  all  the  sections 
SECTIONS  TO  PLOT  -  Each  number  is  taken  as  a 
section  to  be  plotted.  A  letter  (eg.  "S")  will 
end  the  list. 
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If  G,  8,  A  or  a  number  have  not  been  entered  as 
the  first  Input,  control  is  transferred  to  4  to 
try  again. 


VELS 

‘■NO.  OF  FILTERS,  F  -  START,  F  -  END,  Q's" 

B  -  Go  back  to  4. 

NO.  OF  FILTERS  -  The  number  of  filters  to  be  used 
between  frequency  limits  specified  below. 

F  START  -  Starting  frequency. 

F  ENO  -  End  frequency. 

Q  -  Q  of  each  filter. 

The  interval  (F  START  -  F  END)  is  divided  into  a 
number  of  equally  spaced  points  (the  numberof 
filters)  and  a  filter  of  width  F  WIDTH  centered 
at  each  point. 


VELS 

"SYSTEM  FUNCTION  FILE  (N  FOR  NONE)  " 

This  is  the  instrument  section. 

B  -  Go  back  to  5. 

N  -  Do  not  add  instrument 

FILE  NAME  -  File  name  of  instrument  data.  The 
form  is  shown  below; 
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A  REALIO  file  of  the  form: 


Pj^  «  n  the  degrees  of  the  numerator  and 
P2  a  m  denominator 

*^3  *  ‘’n+3  ^o‘  ^1*  *  ®n 

'^n+A  *  ^0’  ^1’  ’  ^m 

where 

ap.  +  a,x^  +  a-x^  + _  +  a  x'^ 

POLYNOMIAL  -  ^ - ^^r-- 

bo  +  b^x  +  b^x  +  _ +  b|^x 

GRPV 

■'SEC  #  ,  OPTION  (H  «  HELP)  " 

If  the  first  input  is  the  letter: 

B  -  Go  back  to  6. 

G  -  Go  on  without  narrow  band  filtering, 
plotting,  etc. 

H  -  List  all  the  available  options. 

S  -  Set  ACUT  and  TCUT  (see  8  for  details). 

F  -  Find  group  velocity  curves  for  the  last 

section  processed. 

C  -  For  each  section,  do  the  narrow  band 

filtering  and  find  the  peaks.  Then  find  the 
group  velocity  curves. 


9 


9 


9 


> 


If  the  first  input  is  the  section  number,  the  options  are: 

N  -  Narrow  band  filter,  etc.  then  plot  the  group 
velocity  curve. 

P  -  Do  a  group  velocity  plot  on  the  data  as  is. 

L  ~  List  the  velocities  for  this  section  number. 

E  -  (Edit)  Load  in  new  values  for  the  frequencies 
and  group  velocities. 

M  -  (Modify)  Use  cross  hairs  to  change  group 
velocities. 


8  GRPV 

“AMP  CUT,  TIM  CUT  ” 

This  request  is  made  if  "S”  has  been  specified  in 

7. 

AMP  CUT  -  Amplitude  cutoff  for  the  group 
velocity  finder  section. 
time  cut  -  Time  cutoff  for  the  group  velocity 
finder  section. 


> 


X 


9  GRPV 

“INDEX,  FREQ,  GRP  V  “ 

This  prompt  is  made  if  "E"  was  specified  in  7. 

This  option  allows  the  modification  of  existing 
data.  Three  data  must  be  entered  to  complete 
each  change,  after  which  the  edit  made  can  be 
terminated  by  entering  any  letter. 

INDEX  -  The  index  of  the  set  to  be  changed. 

FREQ  -  The  frequency. 

GRP  V  -  The  group  velocity. 

137 


S-CL’BED 


If  "M"  is  specified  in  7,  then  cross  hairs  are 
used.  Specify  index  only.  When  cross  hairs 
appear,  adjust  to  desired  point  and  press  any 
number  to  enter.  End  by  setting  INDEX  »  0. 

VELS 

"ENTER  PHASE  CORRCTN  IN  MULTIPLES  OF  PI  (H  =  HELP)" 

5  -  Go  back  to  7  and  do  the  GRPV  section  over. 

H  -  Print  the  options  as  a  user  aid. 

PHACOR  -  (Phase  correction;  will  be  multiplied  by  ir) 
For  the  case  of  step  explosion,  Rayleigh 
wave,  vertical  component,  up  is  positive  and 
PHACOR  is  entered  as  -  3/4.  This  value  will 
be  multiplied  by  ir. 

PHAV 

"GRP  V  SECT  FOR  MATCHED  FILTER  (A  -  1  ON  1)  " 

B  -  Go  back  to  10 

A  -  For  each  section,  the  group  velocity  for  the 
matched  filter  is  that  of  the  section. 

GRPV  -  Input  group  velocity  for  the  matched  filter. 

PHAV 

"MATCHED  FILTER  8W  -  (W)HITE,  (F)ILTER  RANGE  " 

6  -  Go  back  to  II. 

W  -  The  frequencies  go  from  zero  to  the  Nyquist 
frequency. 

F  -  The  frequencies  go  from  zero  to  the  highest 
filter  frequency. 

BW  -  The  frequencies  go  from  zero  to  BW. 


PHITER  (ENTRY  POINT  IN  PHAV) 

"WINDOW  LIMITS  (N  FOR  NONE)?" 

6  -  Go  back  to  12. 

N  -  No  window  limits. 

WINDOW  LIMITS  -  The  lower  and  upper  window  limits  are  read. 
MPLOT 

"SECTIONS  TO  PLOT  -  S  »  ENDLIST,  A  -  ALL,  G  -  GO" 

MPLOT  plots  the  matched  filter  generated  in  PHAV. 

G  -  Go  on  without  plotting. 

A  -  Plot  all  sections 

SECTIONS  TO  PLOT  -  Read  sections  to  plot  until  an 
"S"  Is  read. 

B  -  Go  back  to  13. 

XPLQT 

"SECTIONS  TO  PLOT  -  S  -  ENDLIST,  A  .  ALL,  G  .  GO" 

XPLQT  plots  the  cross-correlated  matched  filter  outputs. 


G  -  Go  on  without  plotting. 

8  -  Go  back  to  14. 

A  -  Plot  all  sections. 

SECTIONS  TO  PLOT  -  Read  section  numbers  to  plot 
until  an  "S"  is  encountered. 
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16  PHITER  (ENTRY  POINT  IN  PHAV) 

••PHASE  VEL  PLOT  -  SECT#  ,  N-TWOPI  “ 

B  -  Go  back  to  15. 

G  -  Go  on  without  plotting, 
any  letter  -  Go  back  to  prompt  at  16. 

(SECT  4  ,  N-TWOPI)  where: 

SECT  #  -  The  section  number  to  be  plotted. 
N-TWOPI  -  The  number  of  2ir's  to  add  to  the 
phase. 

17  PHITER  (ENTRY  POINT  IN  PHAV) 

"OPTION  (H  FOR  HELP)" 

H  -  Point  out  options  available. 

B  -  Go  back  to  15. 

P  -  Loop  over  -  2it,  Oir,  2it,  plotting  the  0i»  with 
a  solid  line  and  the  others  with  a  dotted 
line.  (PPLOT) 

N  -  Make  a  new  phase  velocity  plot.  (PVPLOT) 

G  -  Go  on  without  plot, 

(SECTION  NUMBER,  NTWOPI)  -«  This  option  requires  a 
pair  of  numbers,  the  section  number  and  the  number 
of  2ir's  to  add  to  the  phase  and  increase  phase 
velocity. 

L  -  List  the  frequencies  and  phase  velocities. 


18  GPL0T2 

"SECTION  #,  OPTION  (G  -  GO,  N  .  NEWPLOT)?" 

GPL0T2  plots  the  old  and  new  group  velocity 
curves  to  see  if  a  new  iteration  is  necessary. 
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G  '  Go  on  without  plotting. 

N  -  Go  back  to  the  start  of  the  subroutine  and 
start  a  new  plot. 

B  -  Go  back  to  original  call  to  PHAV  or  PHITER 
(depending  on  the  iteration  number).  The 
first  prompt  following  this  transfer  is  12. 

L  -  List  the  following  data: 

index,  frequency,  old  velocity,  new  velocity. 
At  this  point,  if  the  response  to  the  original 
prompt  was  any  letter,  control  is  transferred 
back  to  the  prompt  at  18. 

SECTION  -  The  section  to  plot. 

At  the  end  of  GPL0T2,  control  is  transferred  back 
to  the  prompt  at  18. 


19  GVITER 

"WANT  NEW  ITERATION  (Y/N)  " 

B  -  Go  back  just  prior  to  19  and  plot  the  old  and 
new  group  velocity  curves  to  see  if  another 
iteration  is  needed. 

N  -  Go  to  prompt  20  and  set  flag  so  new  iteration 
won't  be  done. 

Y  -  Go  to  13,  flag  remains  set  for  new  iteration. 

20  GVITER 

"WANT  THE  NEW  OR  OLD  GROUP  VELOCITES  (N/O)  " 

B  >  Go  back  to  19. 
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N  -  Use  the  new  group  velocites. 
0  -  Use  the  old  group  velocites. 
anything  else  -  Go  back  to  20. 


21  POUT 

“CLEAR  SCREEN  FOR  RUN  SUMMARY" 

This  is  merely  a  command  to  clear  the  screen.  No 
other  action  is  meaningful. 


22  POUT 

“SECTION  TO  PRINT,  N  -  TWOPI  (G  -  GO)?" 

8  -  Go  back  to  call  to  GVITER.  Go  back  to  the 
original  call  to  PHAV  or  PHITER  (depending  on 
the  iteration  number).  The  first  prompt 
following  this  transfer  is  12. 

G  -  Go  on  without  printing, 
any  other  letter  -  Go  to  22  and  redo. 

(SECTION  TO  PRINT,  N-TWO  PI)  -  Specify  the 
section  to  point  and  the  number  of  2it's  to  add  to 
the  phase. 


23  POUT 

“START  FREQ,  END  FREQ,  AND  FREQUENCY  INTERVAL?" 

B  -  Go  back  to  22. 

(SF,  EF,  OF)  where 
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SF  »  Start  frequency. 

EF  a  End  frequency. 

OF  a  Frequency  interval. 


If  any  letter  otner  than  "8“  has  been  entered,  go 
to  23. 


24  POUT 

"OUTPUT  FILE  (N  a  NONE)?" 

N  -  Go  back  to  22. 
3  -  Go  back  to  22. 
Output  file  name. 


25  POUT 

"HEADER  LINE?" 

Anything  entered  (up  to  80  characters)  is  used  as 
a  header  on  the  output  file. 


26  POUT 

"ENTER  L(LOVE)  OR  R(RAYLEIGH)" 

Either  an  L  or  R  is  written  on  the  out  put  file 
to  identify  the  type  of  problem. 

Execution  for  the  given  seismogram  is  complete.  Go  back  to  1  and 
decide  whether  to  start  over  with  a  new  seismogram  or  quit. 
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Figure  1.  Block  diagram,  dispersion  curve  extraction.  The  program  TELVEL  uses  phase-matched  filtering 
as  proposed  by  Herrin  and  Goforth  (1977). 
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